3. Equilibration
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 stage.
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 temperature.
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:
&cntrl
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,
ntb=1, ntp=0,
ntc=2, ntf=2,
nrespa=2,
&end
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).
Line 2:
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.
Line 3:
we will use periodic boundary with constant volume (ntb=1), and no pressure
control (ntp=0)
Line 4:
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
Line 5:
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!
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.