(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 2004
TUTORIAL 7 - SECTION 4
The Nudged Elastic Band Approach to Finding the Lowest
Energy
Pathway Between two States
By Ross Walker
4. Heating the System
The first stage of our NEB MD simulation will be to carefully heat the system up. Note, if you tried to run normal MD with the neb.prmtop and neb.inpcrd files your system would blow up quite spectacularly since these files contain two sets of 15 structures all on top of each other. The electrostatic and VDW repulsion would be infinite on the first step. In the NEB method though the individual molecules never see each other. They only feel a force due to the two springs connecting them to their nearest neightbours.
We are now ready to run MD on our system. The first thing we should do is allow our system to warm up from 0 to 300K. We need to be very careful here though as our initial spring forces will be very large. Hence we do this stage with a small spring constant. Here is the input file we will use, we will run 20ps of MD with a 0.5fs time step (nstlim=40000, dt=0.0005). The half femtosecond time step is probably overkill but it will improve the chances of our simulation remaining stable early on. We will not be using shake to constrain hydrogens in this simulation since the interesting region of the pathway will involve changes in hydrogen bond lengths (ntc=1, ntf=1). This means that even once we have our system equilibrated we will be restricted to a time step of at a maximum 1 fs. Since the coordinates of the end points are fixed in Cartesian space in NEB and the intervening structures cannot move far due to the springs there is no need to recenter the coordinates every few hundred steps so we turn this off (NSCM=0). We will use generalised born for the implicit solvent with a sodium chloride salt concentration of 0.2M (igb=1, saltcon=0.2). The system is small enough here that I will avoid any issues with nonbonded cutoffs and truncated Born radii calculations by setting these to 999.0 angstroms (cut=999.0, rgbmax=999.0). This has the effect of turning off any nonbond cutoff. The NEB specific options are on the last line (ineb=1 [turn on neb], skmin=10, skmax=10 [spring constants]. We start here with relatively small spring constants of 10 KCal/mol/A2. We will increase this value later once the system has begun to equilibrate.
For temperature regulation we use the langevin thermostat (ntt=3) and we will use a very high collision rate of 1000 ps-1 (gamma_ln=1000). This would be extremely high for a regular MD simulation where values of around 5.0 are more appropriate. However, in NEB we are not interested in the actual dynamics, we simply want to ensure that our system explores as much phase space as possible before we cool it back down to 0 K and so we use a high collision rate to ensure rapid uniform heating, give the system some viscosity and improve the stability. Without the strong thermostat, the calculation would become unstable because the periodic motions involving the springs (and the truncation of the forces according to the path tangent) would blow the system apart. Thus it is essential to use very very large values of gamma_ln with NEB simulations.
The large value of gamma_ln would normally mean that our system would heat up to the value of temp0 very rapidly, within a few pico-seconds. This could cause instabilities in the simulation and so we will use NMR weight restraints (nmropt=1) to linearly increase the value of temp0 from 0.0 to 300.0 K over 35,000 steps of the 40000 step simulation (istep1=0,istep2=35000, value1=0.0, value2=300.0). This means that at step 0 the value of the target temperature (temp0) will be 0.0 K. At step 5,000 it will be 42.86 K, at 10,000 it will be 85.72 K, by 20,000 it will be 171.44 K and by 35,000 it will be at 300.0 K. It will then remain at 300.0 K until the end of our 40,000 step run. In this way our system should heat up in a controlled fashion.
|
|
Alanine NEB initial MD with small K &cntrl imin = 0, irest = 0, ntc=1, ntf=1, ntpr=50, ntwx=500, ntb = 0, cut = 999.0, rgbmax=999.0, igb = 1, saltcon=0.2, nstlim = 40000, nscm= 0, dt = 0.0005, ntt = 3, gamma_ln=1000.0, tempi=0, temp0=300, ineb = 1,skmin = 10,skmax = 10, nmropt=1 / &wt type='TEMP0', istep1=0,istep2=35000, value1=0.0, value2=300.0 / &wt type='END' / |
We can now run this MD using the AMBER v9 Sander.PIMD executable (sander.PIMD.MPI if we want to run the simulation in parallel). This will take about 4.5 minutes to run on a single processor of a Pentium-D 3.2GHz Machine:
$AMBERHOME/exe/sander.PIMD -O -i md-ala-neb-smallk_1.in -o md-ala-neb-smallk_1.out -p neb.prmtop -c neb.inpcrd -r md-ala-neb-smallk_1.rst -x md-ala-neb-smallk_1.mdcrd
Input Files: md-ala-neb-smallk_1.in,
neb.prmtop,
neb.inpcrd
Output Files: md-ala-neb-smallk_1.out,
md-ala-neb-smallk_1.rst,
md-ala-neb-smallk_1.mdcrd
Take a look at the output file. It looks just like a normal MD except that you now have the energy of each of the 30 replicates listed. This is the energy of each one as it moves along the potential energy surface. You also see the total energy of the replicates. The lowest energy pathway should be the one for which this sum is the smallest. Thus if we perform simulated annealing calculations and obtain several pathways we can compare the overall sum from each pathway. Take a quick look at the output file and check that the simulation looks stable, i.e. the temperature changes as we expect, slowly increasing to 300.0 K and then oscillating around this value. Any huge temperature jumps are signs of problems with our simulation, poor starting structure, poor choice of simulation options etc. You should also take a look at the replicate energies. These show the potential energy of each replica on the pathway. The most important thing to check at this stage is that on step 0 the first 15 have the same energy, that of the beginning of the pathway and the last 15 have a different but identical energy, that of the end of the pathway. If this is not the case then something is wrong with the neb.prmtop and neb.inpcrd files.
NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 0.00 PRESS = 0.0 Etot = -972.8504 EKtot = 0.0000 EPtot = -992.5236 BOND = 17.6624 ANGLE = 54.1617 DIHED = 210.8264 1-4 NB = 88.0044 1-4 EEL = 1326.1645 VDWAALS = -49.1867 EELEC = -2086.7405 EGB = -553.4159 RESTRAINT = 0.0000 NEB replicate breakdown: Energy for replicate 1 = -33.7654 Energy for replicate 2 = -33.7654 Energy for replicate 3 = -33.7654 Energy for replicate 4 = -33.7654 Energy for replicate 5 = -33.7654 Energy for replicate 6 = -33.7654 Energy for replicate 7 = -33.7654 Energy for replicate 8 = -33.7654 Energy for replicate 9 = -33.7654 Energy for replicate 10 = -33.7654 Energy for replicate 11 = -33.7654 Energy for replicate 12 = -33.7654 Energy for replicate 13 = -33.7654 Energy for replicate 14 = -33.7654 Energy for replicate 15 = -33.7654 Energy for replicate 16 = -32.4029 Energy for replicate 17 = -32.4029 Energy for replicate 18 = -32.4029 Energy for replicate 19 = -32.4029 Energy for replicate 20 = -32.4029 Energy for replicate 21 = -32.4029 Energy for replicate 22 = -32.4029 Energy for replicate 23 = -32.4029 Energy for replicate 24 = -32.4029 Energy for replicate 25 = -32.4029 Energy for replicate 26 = -32.4029 Energy for replicate 27 = -32.4029 Energy for replicate 28 = -32.4029 Energy for replicate 29 = -32.4029 Energy for replicate 30 = -32.4029 Total Energy of replicates = -992.5236 NEB RMS = 1.547934 ------------------------------------------------------------------------------ NMR restraints: Bond = 0.000 Angle = 0.000 Torsion = 0.000 =============================================================================== NSTEP = 50 TIME(PS) = 0.025 TEMP(K) = 0.98 PRESS = 0.0 Etot = -979.3091 EKtot = 1.9186 EPtot = -981.2277 BOND = 24.1932 ANGLE = 58.2929 DIHED = 212.3182 1-4 NB = 88.2867 1-4 EEL = 1324.0833 VDWAALS = -49.3134 EELEC = -2087.5767 EGB = -551.5120 RESTRAINT = 0.0000 NEB replicate breakdown: Energy for replicate 1 = -33.7654 Energy for replicate 2 = -33.7510 Energy for replicate 3 = -33.7567 Energy for replicate 4 = -33.7551 Energy for replicate 5 = -33.7535 Energy for replicate 6 = -33.7561 Energy for replicate 7 = -33.7562 Energy for replicate 8 = -33.7512 Energy for replicate 9 = -33.7539 Energy for replicate 10 = -33.7550 Energy for replicate 11 = -33.7589 Energy for replicate 12 = -33.7547 Energy for replicate 13 = -33.7536 Energy for replicate 14 = -33.7530 Energy for replicate 15 = -27.9942 Energy for replicate 16 = -27.1873 Energy for replicate 17 = -32.3828 Energy for replicate 18 = -32.3861 Energy for replicate 19 = -32.3918 Energy for replicate 20 = -32.3893 Energy for replicate 21 = -32.3946 Energy for replicate 22 = -32.3918 Energy for replicate 23 = -32.3902 Energy for replicate 24 = -32.3863 Energy for replicate 25 = -32.3904 Energy for replicate 26 = -32.3895 Energy for replicate 27 = -32.3935 Energy for replicate 28 = -32.3920 Energy for replicate 29 = -32.3908 Energy for replicate 30 = -32.4029 Total Energy of replicates = -981.2277 NEB RMS = 3.078592 ------------------------------------------------------------------------------ |
We have to do quite a bit of simulated annealing, however, to find the lowest pathway. We also need to increase the strength of our initial 'weak' spring constant.
Before we move on though it is interesting to look at the mdcrd file in vmd. Here you will see all of the structures. Run VMD and load the neb.prmtop and the md-ala-neb-smallk_1.mdcrd file (note: I gzipped the mdcrd file to save space. Ensure you unzip it before loading it in VMD)
Notice how all the structures are still grouped around one of the starting structures. This is because our spring constant wasn't large enough to pull them over the hill between the two end points. It has removed some of the strain on the system however. We are now ready to move onto the next stage of the simulation where we will conduct the simulated annealing and equilibration.