Amber Workshop - tutorial 8

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2005

TUTORIAL 8 - SECTION 5

Case Study: All Atom Structure Prediction and Folding Simulations
of a Stable Protein (Folding Trp-Cage Peptide)

By Ross Walker

Stage 5: Production MD

The final stage of the simulation part of this tutorial is to run a long production simulation at 325 K. In the paper they run 50 ns of simulation at 325 K, but present results for only the first 20 ns, simply stating that the remaining 30 ns showed no significant change from the first part of the simulation. We shall run a similar 50 ns simulation even though the authors saw that the system folded within 5 ns to 20 ns:

"Two independent simulations converged to essentially identical families of structures after 5 and 20 ns."

We will run this 50 ns simulation in 10 stages of 5 ns each. The reason we run it in 10 separate stages is this makes it much easier to resume the simulation should we experience a system crash. It also keeps the output and mdcrd files to manageable sizes. We shall use an identical input file for all 10 stages. The input file is as follows:

equil.in

Stage 2 equilibration 1 0-5ns
 &cntrl
  imin=0, irest=1, ntx=5,
  nstlim=2500000, dt=0.002,
  ntc=2, ntf=2,
  ntt=1, tautp=0.5,
  tempi=325.0, temp0=325.0,
  ntpr=500, ntwx=500,
  ntb=0, igb=1,
  cut=999.,rgbmax=999.
 /

Each of the 10 stages will consist of 2,500,000 steps (nstlim) with a 2 fs time step (dt) giving 5 ns per stage. Note we have SHAKE on the whole time (ntc=2, ntf=2), we use the Berendsen thermostat for temperature control (ntt=1) and now that we have heated our system up and it appears to be stable we can use a much more closely coupled thermostat of 0.5 ps (tautp=0.5). This will serve to keep our system close to 325 K. We set our target temperature at 325K as they do in the paper and we write to the output file and mdcrd file every 500 steps. Writing more frequently that this would result in huge files. The mdcrd file for each 5 ns is 35 Mb as it is. Every 500 steps should be frequent enough for our purposes.

Here is the PBS script I used to run the 10 jobs, you can modify it for your system:

#PBS -l ncpus=16
#PBS -l walltime=500:00:00
#PBS -l cput=8000:00:00
#PBS -j oe
setenv AMBERHOME /usr/people/rcw/amber9
cd ~rcw/production
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c heat7.rst -r equil1.rst -o equil1.out -x equil1.mdcrd
gzip -9 equil1.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil1.rst -r equil2.rst -o equil2.out -x equil2.mdcrd
gzip -9 equil2.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil2.rst -r equil3.rst -o equil3.out -x equil3.mdcrd
gzip -9 equil3.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil3.rst -r equil4.rst -o equil4.out -x equil4.mdcrd
gzip -9 equil4.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil4.rst -r equil5.rst -o equil5.out -x equil5.mdcrd
gzip -9 equil5.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil5.rst -r equil6.rst -o equil6.out -x equil6.mdcrd
gzip -9 equil6.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil6.rst -r equil7.rst -o equil7.out -x equil7.mdcrd
gzip -9 equil7.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil7.rst -r equil8.rst -o equil8.out -x equil8.mdcrd
gzip -9 equil8.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil8.rst -r equil9.rst -o equil9.out -x equil9.mdcrd
gzip -9 equil9.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil9.rst -r equil10.rst -o equil10.out -x equil10.mdcrd
gzip -9 equil10.mdcrd
echo "DONE"

All 10 stages take a total of about 27 hours to run on 16 cpus of a 1.3GHz SGI Altix. Here are the output files, you can download them individually or as a single tar.gz file (Warning the mdcrd files are over 13Mb each compressed):

Heating Stage	Output File	Restrt File	Mdcrd File
Stage 1	equil1.out	equil1.rst	equil1.mdcrd.gz
Stage 2	equil2.out	equil2.rst	equil2.mdcrd.gz
Stage 3	equil3.out	equil3.rst	equil3.mdcrd.gz
Stage 4	equil4.out	equil4.rst	equil4.mdcrd.gz
Stage 5	equil5.out	equil5.rst	equil5.mdcrd.gz
Stage 6	equil6.out	equil6.rst	equil6.mdcrd.gz
Stage 7	equil7.out	equil7.rst	equil7.mdcrd.gz
Stage 8	equil8.out	equil8.rst	equil8.mdcrd.gz
Stage 9	equil9.out	equil9.rst	equil9.mdcrd.gz
Stage 10	equil10.out	equil10.rst	equil10.mdcrd.gz
Complete file set	production.tar.gz (132 Mb)

You should load the mdcrd files into VMD (unzip them first, it may take a while to load them. Perhaps one day the VMD authors will add support for reading compressed files <hint hint>). You can then watch the complete trajectory. It is best if you choose ribbons as the representation. Here is one such structure coloured by residue type:

CLICK HERE TO GO TO SECTION 6