(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 4

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

By Ross Walker

Stage 4: Heating the system up.

The next stage is to do some MD on this system. We will first by slowly heat the system up over 50 ps in a total of 7 stages. Doing the heating in stages reduces the chances that the system will blow up by allowing it to equilibrate at each temperature. An alternative option would be to use weight restraints to control the heating. See the AMBER manual for more information. Normally we would run MD simulations at 300 K (room temperature). However, we should refer to a comment they make in the paper:

MD simulations of 100 ns were performed at 300 K, but all were kinetically trapped on this time scale, showing strong dependence on initial conditions and failing to converge to similar conformational ensembles. We therefore increased the temperature to 325 K.

They found that it was necessary to heat the system to 325 K to avoid being kinetically trapped in local minima. We shall therefore do the same. When running at elevated temperatures you need to think carefully about the timestep since oscillations in the system will be more pronounced. This means that if you were to run at 600 K with a time step of 2 fs, that would normally be sufficient for a shaken system at 300 K, the time between updates would be too long and lead to instabilities. Fortunately 325 K is sufficiently close to 300 K that we can get away with a 2 fs time step as long as we constrain bonds involving hydrogen. However, if we later heat the system to 400 K, we would likely need to reduce the time step to 1.5 fs or so.

Since our initial structure is a manually built structure rather than an experimental crystal structure, it is likely to initially be much less stable than an experimental structure. To allow our initial system to relax in a controlled fashion, we will use a very short time step of 0.5 fs for the heating stage and then switch to 2 fs for the equilibration. This is probably overkill but it is better to be safe than sorry. We will also set the coordinate (mdcrd) write frequency to 50 here (ntwx = 50). This is a very high write frequency and could in large parallel simulations hurt performance but during the very initial heating there is an increased chance of system instability, especially since this is a model built structure. Setting a low number of steps for each coordinate write means that if problems do occur we should be able to see what happens and identify where the problem comes from. When it comes to doing production MD more appropriate values are in the range of every 500 to 1000 steps. This system has only 304 atoms and so each coordinate frame will be quite small so we can safely write every 500 steps or so in the production phase without impacting performance. However, if this was a large explicit solvent simulation running in parallel on a large cluster then, unless more frequent sampling is required for some reason, we should aim to write the coordinates every 2 to 5 ps or so. For a 2fs time step this would mean setting ntwx between 1000 and 2500.

So, our protocol for heating the system will be as follows:

    Step 1 - 10,000 steps, 0.5fs time step (5ps), initial temp = 0.0K, target temp = 50.0K, temperature coupling constant = 1.0ps
    Step 2 - 10,000 steps, 0.5fs time step (5ps), target temp = 100.0K, temperature coupling constant = 1.0ps
    Step 3 - 10,000 steps, 0.5fs time step (5ps), target temp = 150.0K, temperature coupling constant = 1.0ps
    Step 4 - 10,000 steps, 0.5fs time step (5ps), target temp = 200.0K, temperature coupling constant = 1.0ps
    Step 5 - 10,000 steps, 0.5fs time step (5ps), target temp = 250.0K, temperature coupling constant = 1.0ps
    Step 6 - 10,000 steps, 0.5fs time step (5ps), target temp = 300.0K, temperature coupling constant = 1.0ps
    Step 7 - 40,000 steps, 0.5fs time step (20ps), target temp = 325.0K, temperature coupling constant = 1.0ps

Here is the first input file:

Stage 1 heating of TC5b 0 to 50K
 &cntrl
  imin=0, irest=0, ntx=1,
  nstlim=10000, dt=0.0005,
  ntc=2, ntf=2,
  ntt=1, tautp=1.0,
  tempi=0.0, temp0=50.0,
  ntpr=50, ntwx=50,
  ntb=0, igb=1,
  cut=999.,rgbmax=999.
 /

The other 6 input files are very similar, with just the relevant parameters changed. They are available here: (heat2.in, heat3.in, heat4.in, heat5.in, heat6.in, heat7.in).

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

#PBS -l ncpus=16
#PBS -l walltime=500:00:00
#PBS -l cput=2000:00:00
#PBS -j oe
setenv AMBERHOME /usr/people/rcw/amber9
cd ~rcw/initial_heating
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrd
gzip -9 heat1.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat2.in -p TC5b.prmtop -c heat1.rst -r heat2.rst -o heat2.out -x heat2.mdcrd
gzip -9 heat2.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat3.in -p TC5b.prmtop -c heat2.rst -r heat3.rst -o heat3.out -x heat3.mdcrd
gzip -9 heat3.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat4.in -p TC5b.prmtop -c heat3.rst -r heat4.rst -o heat4.out -x heat4.mdcrd
gzip -9 heat4.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat5.in -p TC5b.prmtop -c heat4.rst -r heat5.rst -o heat5.out -x heat5.mdcrd
gzip -9 heat5.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat6.in -p TC5b.prmtop -c heat5.rst -r heat6.rst -o heat6.out -x heat6.mdcrd
gzip -9 heat6.mdcrd
mpirun -np 16 $AMBERHOME/exe/sander -O -i heat7.in -p TC5b.prmtop -c heat6.rst -r heat7.rst -o heat7.out -x heat7.mdcrd
gzip -9 heat7.mdcrd
echo "DONE"

These 7 simulations take a total of around 7 mins on 16 processors of a 1.3GHz SGI Altix. Here are the output files, you can download them individually or as a single tar.gz file:

You should load the mdcrd files into VMD (unzip them first) so that you can watch what happens during the heating. You should see the system start to fold up. We are not too worried about the structures that are sampled at this stage, we just need to verify that things look okay and nothing strange is happening.

Here is the final structure of heat7.mdcrd:

We are beginning to see some alpha helix formation but we have a long way to go before we start seeing stable folded structures.

CLICK HERE TO GO TO SECTION 5

(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