Amber masthead
Filler image AmberTools24 Amber24 Manuals Tutorials Force Fields Contacts History
Filler image

Useful links:

Amber Home
Download Amber
Installation
Amber Citations
GPU Support
Updates
Mailing Lists
For Educators
File Formats
Contributors
Workshops
Sampling Configuration Space
 

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

Formerly known as TUTORIAL A5: The Nudged Elastic Band Approach to Finding the Lowest Energy Pathway Between two States

By Christina Bergonzo, Carlos Simmerling & Ross Walker

Return to Index Page

6. Slow Cooling

The final stage is to slowly cool the system down so that the images freeze out along the pathway. We will do the cooling in 2 stages. In the first stage we shall gradually cool down in 50 K steps over 120 ps using the following profile:

0 to 10ps : 300 K -> 250 K
10 to 20ps : Equil at 250 K
20 to 30ps : 250 K -> 200 K
30 to 40ps : Equil at 200 K
40 to 50ps : 200 -> 150 K
50 to 60ps : Equil at 150 K
60 to 70ps : 150 K -> 100 K
70 to 80ps : Equil at 100 K
80 to 90ps : 100 K -> 50 K
90 to 100ps : Equil at 50 K
100 to 110 ps : 50 K -> 0 K
110 to 120 ps : Equil at 0 K

Note since we are now back at 300 K and below we can safely switch back to using a 1 fs time step. Here is the script to create the input file and groupfile:

groupgen-cool.sh

#!/bin/bash

mkdir 4COOL1/
cd 4COOL1/

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

cat > 4cool1.in <<EOF
120ps NEB ALA-ALA cool stage 1
 &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 = 120000, 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",
  nmropt=1,ig=-1,
 /
 &wt type='TEMP0', istep1=0,istep2=10000,
   value1=300.0, value2=250.0
 /
 &wt type='TEMP0', istep1=10001,istep2=20000,
   value1=250.0, value2=250.0
 /
 &wt type='TEMP0', istep1=20001,istep2=30000,
   value1=250.0, value2=200.0
 /
 &wt type='TEMP0', istep1=30001,istep2=40000,
   value1=200.0, value2=200.0
 /
 &wt type='TEMP0', istep1=40001,istep2=50000,
   value1=200.0, value2=150.0
 /
 &wt type='TEMP0', istep1=50001,istep2=60000,
   value1=150.0, value2=150.0
 /
 &wt type='TEMP0', istep1=60001,istep2=70000,
   value1=150.0, value2=100.0
 /
 &wt type='TEMP0', istep1=70001,istep2=80000,
   value1=100.0, value2=100.0
 /
 &wt type='TEMP0', istep1=80001,istep2=90000,
   value1=100.0, value2=50.0
 /
 &wt type='TEMP0', istep1=90001,istep2=100000,
   value1=50.0, value2=50.0
 /
 &wt type='TEMP0', istep1=100001,istep2=110000,
   value1=50.0, value2=0.0
 /
 &wt type='TEMP0', istep1=110001,istep2=120000,
   value1=0.0, value2=0.0
 /
 &wt type='END'
 /
EOF
cd ..

The file creates a directory named 4COOL1, since this is the first cooling step. We can cd to this directory and run sander:

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

This takes about 22 minutes, and produces the following output files: 4COOL1.tar.gz

Normally we would be done running simulated annealing at this point, however, the NEB method needs a final energy minimization of the path. This is done here using the velocity Verlet algorithm:

 

We will do one final 200ps long stage of cooling where we will set temp0=0.0 K and we will turn on the quenched velocity verlet method that was discussed in section 1 (vv=1). Here is the last script which creates an input file and groupfile for the final quenched MD at 0K:

groupgen-cool2.sh

#!/bin/bash

mkdir 5COOL2/
cd 5COOL2/

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

cat > 5cool2.in <<EOF
50ps NEB ALA-ALA cool stage 2
 &cntrl
  imin = 0, irest = 1, ntx=5,
  ntc=1, ntf=1,
  ntpr=100, ntwx=10000,
  ntb = 0, cut = 999.0, rgbmax=999.0,
  igb = 1, saltcon=0.2,
  nstlim = 200000, nscm= 0,
  dt = 0.001,
  ntt = 3, gamma_ln=1000.0,
  temp0=0.0,
  ineb = 1,skmin = 50,skmax = 50,
  tgtfitmask=":1,2,3",
  tgtrmsmask=":1,2,3@N,CA,C",
  vv=1,vfac=0.1,ig=-1,
 /
EOF
cd ..

We can run this as we do above:

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

This takes about 40 minutes, and produces the following output files: 5COOL2.tar.gz

Lets take a look at the final step of this final cooling (looking at neb.out.003):

   NSTEP =   200000   TIME(PS) =     740.000  TEMP(K) =     0.00  PRESS =     0.0
 Etot   =       -36.1351  EKtot   =         0.0000  EPtot      =       -36.1351
 BOND   =         0.5727  ANGLE   =         1.3399  DIHED      =         8.9864
 1-4 NB =         2.9599  1-4 EEL =        44.7642  VDWAALS    =        -1.8926
 EELEC  =       -75.0851  EGB     =       -17.7805  RESTRAINT  =         0.0000
NEB replicate breakdown:
Energy for replicate   1 =      -36.0518
Energy for replicate   2 =      -36.1171
Energy for replicate   3 =      -36.1171
Energy for replicate   4 =      -36.1351
Energy for replicate   5 =      -36.1475
Energy for replicate   6 =      -36.1528
Energy for replicate   7 =      -36.1528
Energy for replicate   8 =      -36.1529
Energy for replicate   9 =      -36.1428
Energy for replicate  10 =      -36.1130
Energy for replicate  11 =      -36.1129
Energy for replicate  12 =      -36.1129
Energy for replicate  13 =      -36.0827
Energy for replicate  14 =      -36.0490
Energy for replicate  15 =      -36.0064
Energy for replicate  16 =      -35.9555
Energy for replicate  17 =      -35.8801
Energy for replicate  18 =      -35.8376
Energy for replicate  19 =      -35.6213
Energy for replicate  20 =      -35.4317
Energy for replicate  21 =      -34.9259
Energy for replicate  22 =      -33.9905
Energy for replicate  23 =      -32.3218
Energy for replicate  24 =      -30.4850
Energy for replicate  25 =      -29.3411
Energy for replicate  26 =      -29.8267
Energy for replicate  27 =      -31.2547
Energy for replicate  28 =      -32.0472
Energy for replicate  29 =      -32.7606
Energy for replicate  30 =      -33.5716
Energy for replicate  31 =      -33.9966
Energy for replicate  32 =      -34.1210
Total Energy of replicates =    -1109.0158

 

 

Notice how the energies fluctuate along the pathway. It is probably easier to visualize this if we plot the energies. We can manually extract them an place them in a file: final_pathway_energies.dat.

As you can see there is 1 major maxima along this pathway and the overall barrier height is approximately 7 KCal/mol in one direction and 4 KCal/mol in the reverse direction. In reality you would probably want to use at least 32 images here to get a nice smooth curve.