In MD simulations, atoms of the macromolecules and of the surrounding solvent
undergo a relaxation that usually lasts for tens or hundreds of picoseconds
before the system reaches a stationary state. The initial nonstationary
segment of the simulated trajectory are typically discarded in the calculation
of equilibrium properties. This stage of the MD simulation is called
Equilibration protocols are still largely a matter of personal preference.
Some protocols call for very elaborate procedures involving gradually increasing
temperature in a step-wise fashion while other more aggressive approach simply
use a linear temperature gradient and heat the system up to the desired
In our example, we'll follow the protocol suggested in the AMBER manual and
perform a two stage equilibration. In the first stage, we will start the
system from a low temperature of 100 K and gradually heat up to 300 K over 10
picosecond of simulation time. We will perform this stage of equilibration
with the volume held constant. The input file is as follows:
Heating up the system equilibration stage 1
nstlim=5000, dt=0.002, ntx=1, irest=0, ntpr=500, ntwr=5000, ntwx=5000,
tempi =100.0, temp0=300.0, ntt=1, tautp=2.0, ig=209858,
Let's go through this line by line.
we specify that we want to run for 5000 steps (nstlim=5000), with a time step
of 2 fs (dt=0.002); we are starting a new run so don't have previous velocity
information (ntx=1, irest=0), we will print energy output every 5000 steps (ntwr=5000)
and save coordinates every 5000 (ntwx=5000).
we want the initial temperature to be 100K (tempi=100.0) and the reference
temperature to be 300K (temp0=300.0) using Berendsen coupling algorithm to
maintain constant temperature. The time constant for temperature
coupling is 2 ps (tautp=2.0). For the initial velocity assignment, a
random number seed is required (ig=209858). If you use the same seed,
you will get the exact same trajectory. If you want to run multiple
trajectories with different assignment of initial velocities, you need to use
a different random seed number.
we will use periodic boundary with constant volume (ntb=1), and no pressure
we will use the SHAKE algorithm to constrain only bonds involving hydrogen (ntc=2)
and omit force evaluations for these bonds (ntf=2). We also set a
tighter tolerance from resetting coordinates
to speed up the calculations a bit, we also take a bigger time step (nrespa=2)
for evaluating the slow-varying terms in the force field. In this
setting, the time step is equal to nrespa * dt ( 2 * 2 fs = 4 fs.)
note: the nrespa option only works with AMBER 7 and above. If you are using
older versions of AMBER, remember to exclude nrespa in your input.
Ok, save this input file as "eq_v.in". Now let's fire up the engine and
give it a go!
$ sander -O -i eq_v.in -p wt1mg.parm7 -c wt1mg_min_all.rst -r wt1mg_eq_v.rst
-x wt1mg_eq_v.crd -o wt1mg_eq_v.out
This calculation took 2 hours on a dual cpu Pentium III, so you can take
another coffee break here. Alternatively, you can download my results from
the download page.