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

Nudged Elastic Band (NEB) simulations - SECTION 5

By Christina Bergonzo, Carlos Simmerling & Ross Walker

Return to Index Page

5. Simulated Annealing and Equilibration

We will now run 100ps of MD at 300K with a much bigger spring constant. This will allow the structures to explore the energy landscape. At the same time we should be able to increase the time step to 1 fs now that we have successfully heated the system to 300 K and allowed it to relax slightly. Since we are only really going to be interested in the overall pathway, rather than the dynamics as we would be with a regular MD simulation I will set the mdcrd frequency to once every 5,000 steps (ntwx=5000). This allows us to avoid using too much disk space for a trajectory file we won't be using. We could turn of trajectory writing altogether but I like to have it write to the trajectory file occasionally as it can help in locating problems if things go wrong.

Similarly to the previous step, we will use a script to make a directory (here named 2EQUIL), create an input file as described above, and create a groupfile. The script is below, called groupgen-equilibrate.sh.

groupgen-equilibrate.sh

#!/bin/bash
#Note: If you copy this script verbatim, take out the space before each of the "EOF"s before using.

mkdir 2EQUIL/
cd 2EQUIL/

if [[ -e groupfile ]] ; then
  rm groupfile
fi
for ((i=0; i <=31 ; i++)) ; do
 ext=`printf "%03i" $i`
 echo -O -p ../str1.prmtop -c ../1HEAT/neb.r.$ext -i ./2equil.in -x ./neb.x.$ext \
 -o ./neb.out.$ext -inf ./neb.info.$ext -r ./neb.r.$ext >> groupfile

done
cat > 2equil.in << EOF
100ps NEB ALA-ALA equilibration
&cntrl
 imin = 0, irest = 1, ntx=5,
 ntc=1, ntf=1,
 ntpr=1000, ntwx=5000,
 ntb = 0, cut = 999.0, rgbmax=999.0,
 igb = 1, saltcon=0.2,
 nstlim = 100000, nscm= 0,
 dt = 0.001,
 ntt = 3, gamma_ln=1000.0,
 temp0=300,
 ineb = 1,skmin = 50,skmax = 50
 tgtfitmask=":1,2,3",
 tgtrmsmask=":1,2,3@N,CA,C",
 ig=-1,
/

 EOF
 
cd ..

Here is the input file created by the script:

100ps NEB ALA-ALA equilibration
 &cntrl
 imin = 0, irest = 1, ntx=5,
 ntc=1, ntf=1,
 ntpr=1000, ntwx=5000,
 ntb = 0, cut = 999.0, rgbmax=999.0,
 igb = 1, saltcon=0.2,
 nstlim = 100000, nscm= 0,
 dt = 0.001,
 ntt = 3, gamma_ln=1000.0,
 temp0=300,
 ineb = 1,skmin = 50,skmax = 50
 tgtfitmask=":1,2,3",
 tgtrmsmask=":1,2,3@N,CA,C",
 ig=-1,
 /

and here is the groupfile:

-O -p ../str1.prmtop -c ../1HEAT/neb.r.000 -i ./2equil.in -x ./neb.x.000 -o ./neb.out.000 -inf ./neb.info.000 -r ./neb.r.000
-O -p ../str1.prmtop -c ../1HEAT/neb.r.001 -i ./2equil.in -x ./neb.x.001 -o ./neb.out.001 -inf ./neb.info.001 -r ./neb.r.001
-O -p ../str1.prmtop -c ../1HEAT/neb.r.002 -i ./2equil.in -x ./neb.x.002 -o ./neb.out.002 -inf ./neb.info.002 -r ./neb.r.002
-O -p ../str1.prmtop -c ../1HEAT/neb.r.003 -i ./2equil.in -x ./neb.x.003 -o ./neb.out.003 -inf ./neb.info.003 -r ./neb.r.003
-O -p ../str1.prmtop -c ../1HEAT/neb.r.004 -i ./2equil.in -x ./neb.x.004 -o ./neb.out.004 -inf ./neb.info.004 -r ./neb.r.004
-O -p ../str1.prmtop -c ../1HEAT/neb.r.005 -i ./2equil.in -x ./neb.x.005 -o ./neb.out.005 -inf ./neb.info.005 -r ./neb.r.005
-O -p ../str1.prmtop -c ../1HEAT/neb.r.006 -i ./2equil.in -x ./neb.x.006 -o ./neb.out.006 -inf ./neb.info.006 -r ./neb.r.006
-O -p ../str1.prmtop -c ../1HEAT/neb.r.007 -i ./2equil.in -x ./neb.x.007 -o ./neb.out.007 -inf ./neb.info.007 -r ./neb.r.007
-O -p ../str1.prmtop -c ../1HEAT/neb.r.008 -i ./2equil.in -x ./neb.x.008 -o ./neb.out.008 -inf ./neb.info.008 -r ./neb.r.008
-O -p ../str1.prmtop -c ../1HEAT/neb.r.009 -i ./2equil.in -x ./neb.x.009 -o ./neb.out.009 -inf ./neb.info.009 -r ./neb.r.009
-O -p ../str1.prmtop -c ../1HEAT/neb.r.010 -i ./2equil.in -x ./neb.x.010 -o ./neb.out.010 -inf ./neb.info.010 -r ./neb.r.010
-O -p ../str1.prmtop -c ../1HEAT/neb.r.011 -i ./2equil.in -x ./neb.x.011 -o ./neb.out.011 -inf ./neb.info.011 -r ./neb.r.011
-O -p ../str1.prmtop -c ../1HEAT/neb.r.012 -i ./2equil.in -x ./neb.x.012 -o ./neb.out.012 -inf ./neb.info.012 -r ./neb.r.012
-O -p ../str1.prmtop -c ../1HEAT/neb.r.013 -i ./2equil.in -x ./neb.x.013 -o ./neb.out.013 -inf ./neb.info.013 -r ./neb.r.013
-O -p ../str1.prmtop -c ../1HEAT/neb.r.014 -i ./2equil.in -x ./neb.x.014 -o ./neb.out.014 -inf ./neb.info.014 -r ./neb.r.014
-O -p ../str1.prmtop -c ../1HEAT/neb.r.015 -i ./2equil.in -x ./neb.x.015 -o ./neb.out.015 -inf ./neb.info.015 -r ./neb.r.015
-O -p ../str1.prmtop -c ../1HEAT/neb.r.016 -i ./2equil.in -x ./neb.x.016 -o ./neb.out.016 -inf ./neb.info.016 -r ./neb.r.016
-O -p ../str1.prmtop -c ../1HEAT/neb.r.017 -i ./2equil.in -x ./neb.x.017 -o ./neb.out.017 -inf ./neb.info.017 -r ./neb.r.017
-O -p ../str1.prmtop -c ../1HEAT/neb.r.018 -i ./2equil.in -x ./neb.x.018 -o ./neb.out.018 -inf ./neb.info.018 -r ./neb.r.018
-O -p ../str1.prmtop -c ../1HEAT/neb.r.019 -i ./2equil.in -x ./neb.x.019 -o ./neb.out.019 -inf ./neb.info.019 -r ./neb.r.019
-O -p ../str1.prmtop -c ../1HEAT/neb.r.020 -i ./2equil.in -x ./neb.x.020 -o ./neb.out.020 -inf ./neb.info.020 -r ./neb.r.020
-O -p ../str1.prmtop -c ../1HEAT/neb.r.021 -i ./2equil.in -x ./neb.x.021 -o ./neb.out.021 -inf ./neb.info.021 -r ./neb.r.021
-O -p ../str1.prmtop -c ../1HEAT/neb.r.022 -i ./2equil.in -x ./neb.x.022 -o ./neb.out.022 -inf ./neb.info.022 -r ./neb.r.022
-O -p ../str1.prmtop -c ../1HEAT/neb.r.023 -i ./2equil.in -x ./neb.x.023 -o ./neb.out.023 -inf ./neb.info.023 -r ./neb.r.023
-O -p ../str1.prmtop -c ../1HEAT/neb.r.024 -i ./2equil.in -x ./neb.x.024 -o ./neb.out.024 -inf ./neb.info.024 -r ./neb.r.024
-O -p ../str1.prmtop -c ../1HEAT/neb.r.025 -i ./2equil.in -x ./neb.x.025 -o ./neb.out.025 -inf ./neb.info.025 -r ./neb.r.025
-O -p ../str1.prmtop -c ../1HEAT/neb.r.026 -i ./2equil.in -x ./neb.x.026 -o ./neb.out.026 -inf ./neb.info.026 -r ./neb.r.026
-O -p ../str1.prmtop -c ../1HEAT/neb.r.027 -i ./2equil.in -x ./neb.x.027 -o ./neb.out.027 -inf ./neb.info.027 -r ./neb.r.027
-O -p ../str1.prmtop -c ../1HEAT/neb.r.028 -i ./2equil.in -x ./neb.x.028 -o ./neb.out.028 -inf ./neb.info.028 -r ./neb.r.028
-O -p ../str1.prmtop -c ../1HEAT/neb.r.029 -i ./2equil.in -x ./neb.x.029 -o ./neb.out.029 -inf ./neb.info.029 -r ./neb.r.029
-O -p ../str1.prmtop -c ../1HEAT/neb.r.030 -i ./2equil.in -x ./neb.x.030 -o ./neb.out.030 -inf ./neb.info.030 -r ./neb.r.030
-O -p ../str1.prmtop -c ../1HEAT/neb.r.031 -i ./2equil.in -x ./neb.x.031 -o ./neb.out.031 -inf ./neb.info.031 -r ./neb.r.031

You can see in the above groupfile that the starting coordinates are taken from the previous step's final coordinates. This is important, since the most recent coordinates along the path are the best resolved images.

We can now run this stage of MD using the same command as before, except in the 2EQUIL directory:

mpirun -np 32 $AMBERHOME/exe/sander.MPI -ng 32 -groupfile groupfile

This will take about 18 minutes, and will create the following output files: 2EQUIL.tar.gz

The next step is to run a simulated annealing run. Here we shall use the NMR restraints option again to control the value of temp0 during the run. We will run a total of 300ps of MD (this should be sufficient for our small system) with the following temperature profile:

0 - 50ps : Heat from 300 K to 400 K
50 - 100ps : Equilibrate at 400 K
100 - 150ps : Heat from 400 to 500 K
150 - 200ps : Equilibrate at 500 K
200 - 250ps : Cool to 300 K
250 - 300ps : Equilibrate at 300 K

Note, as we will be heating to 500K it is necessary to reduce the time step to 0.5fs. This is because at these higher temperatures a 1fs time step without shake would cause instabilities. This annealing protocol is probably overkill for such a small system but it serves to illustrate how to do it. Again, we will use a script to create the input file and groupfile.

groupgen-anneal.sh

#!/bin/bash
# Note: take out the space before both "EOF" values if you copy this verbatim.

mkdir 3ANNEAL/
cd 3ANNEAL/

if [[ -e groupfile ]] ; then
  rm groupfile
fi

for ((i=0; i <= 31; i++)) ; do
  ext=`printf "%03i" $i`
  echo -O -p ../str1.prmtop -c ../2EQUIL/neb.r.$ext -i ./3anneal.in -x ./neb.x.$ext -o ./neb.out.$ext -inf ./neb.info.$ext -r ./neb.r.$ext >> groupfile
done

cat > 3anneal.in << EOF

300ps NEB ALA-ALA simulated annealing
 &cntrl
  imin = 0, irest = 1, ntx=5,
  ntc=1, ntf=1,
  ntpr=500, ntwx=10000,
  ntb = 0, cut = 999.0, rgbmax=999.0,
  igb = 1, saltcon=0.2,
  nstlim = 600000, nscm= 0,
  dt = 0.0005,
  ntt = 3, gamma_ln=1000.0,
  temp0=300,
  ineb = 1,skmin = 50,skmax = 50,
  tgtfitmask=":1,2,3",
  tgtrmsmask=":1,2,3@N,CA,C",
  nmropt=1, ig=-1,
 /
 &wt type='TEMP0', istep1=0,istep2=100000,
   value1=300.0, value2=400.0
 /
 &wt type='TEMP0', istep1=100001,istep2=200000,
   value1=400.0, value2=400.0
 /
 &wt type='TEMP0', istep1=200001,istep2=300000,
   value1=400.0, value2=500.0
 /
 &wt type='TEMP0', istep1=300001,istep2=400000,
   value1=500.0, value2=500.0
 /
 &wt type='TEMP0', istep1=400001,istep2=500000,
   value1=500.0, value2=300.0
 /
 &wt type='TEMP0', istep1=500001,istep2=600000,
   value1=300.0, value2=300.0
 /
 &wt type='END'
 /

 EOF

cd ..

We can now run this stage of MD using the same command as before, except in the 3ANNEAL directory:

mpirun -np 32 $AMBERHOME/exe/sander.MPI -ng 32 -groupfile groupfile

This takes about 1 hour 45 minutes , and gives the following output files: 3ANNEAL.tar.gz

Before we proceed to the next stage and cool the system down so that it freezes out along the low energy pathway we should quickly check that the temperature profile we got from our simulated annealing run is what we expected. Here I will use grep to extract every line in the output file that contains the word NSTEP. I then 'pipe' this to awk which I use to extract only the 6th and 9th columns from each of those lines (these are the numerical values immediately after TIME = and TEMP =) and I write this to a file called temp_profile.dat. I have deleted the last two lines of this file as these two points refer to the RMS and average values. This should be done for every output file but for the purposes of this tutorial we will just look at the 4th point (neb.out.003).You can then plot this with whichever graphing program you like. As this is just a quick check that things proceeded okay I will use xmgrace:

grep NSTEP neb.out.003 | awk '{print $6 "\t" $9}' > temp_profile.dat

xmgrace temp_profile.dat

As you can see above everything looks okay, the degree of fluctuations at the end of the trajectory, when we have cooled back down to 300 K, are roughly the same as those at the beginning of the trajectory and there are no sudden jumps in the temperature so everything looks good. Now we are ready to progress to the next step.