Re: [AMBER] Different trajectories with different RMSD with changing the parameters of gamma_ln, taup, timesteps and extend of density equilibration

From: Thomas Cheatham III <tec3.utah.edu>
Date: Fri, 8 Jul 2011 15:49:34 -0600 (Mountain Daylight Time)

Re: http://archive.ambermd.org/201107/0123.html
    http://archive.ambermd.org/201107/0093.html

Note that if people do not reply in a timely manner, please do not simply
repost older messages. Try re-phrasing / condensing and/or make it
clear that you are asking again since you did not get sufficient help
previously.

> I am working on DNA-protein interactions. I have followed the
> following steps:
...
> I have performed 9 (A-I) different simulations with different
> temperature coupling parameter (gamma_ln ) and pressure couling
...
> I calculated the RMSD deviations withrest to the first frame (BLACK
> lines in the RMSD plots) in the production stage and also with repect
> to the minimized structure after minimzation step2 (RED lines).. I
> have attached the RMSD plots. The RMSD plots shows a large variation
> upon changing the simulation parameters..?
...
> If these trajectories differ by changing the simulations parameters,
> how we can trust these results?

MD simulations are chaotic and small changes in parameters or initial
conditions (and/or parallelism) will lead to different results. As
pointed out by Professor Roitberg, a key question is whether the observed
differences are expected and/or statistically significant.

Before assuming that the differences observed are real, convince yourself
that they are statistically significant. Roitberg suggested that you
cannot easily compare results from different runtime parameters (i.e.
gamma_ln); you replied that perhaps that meant you should re-run all the
simulations (to which no one replied).

Since you have set ig=-1, re-running should lead to different results
(noting that you may want to set ig=-1 in the heat1 phase). However, you
need not re-run all of them. Pick one set of conditions (gamma_ln, taup,
dt) and run these simulations multiple times. What you likely will
observe is similar differences to what you saw previously under multiple
conditions, i.e. that the system has not converged and the results vary
considerably from simulation to simulation (although I could be wrong
here). Another good way to get different initial conditions is to use a
different ion distribution (you can use randomizeions in ptraj).

> Is there any 'best parameters for gamma_ln and taup.

The less often you perturb the temperature, the better, however I would
not expect significant differences in the results for simulations with
gamma_ln = 1 or 2 or 3. For equilibration, I set taup smaller (0.2 or
1.0) to get the pressure correct quickly, and then during production
either run constant volume or set taup larger (5.0, 10.0) for less
perturbation.

> Is there any problems in the MD protocol which I followed?

The protocol is (mostly) fine; likely the observed differences are not
statistically significant.

> What is the reason for the odd behaviour in the RMSD plots and what
> should I do to get a stable MD trajectory?

A subtlety relates to the restraints on the protein+DNA with constant
pressure since the positions of the molecules are scaled shifting the
relative positions of the molecules as the simulation runs. If the box
size changes significantly (it shrinks by default with standard LEaP
solvateOct/solvateBox), this can lead to the molecules getting closer
together despite the restraints. To get around this, the molecules have
to be merged into a single molecule which requires hacking the prmtop.

  http://archive.ambermd.org/200609/0171.html

If this molecule shifting is leading to instability (i.e. the structural
distortion you are seeing), you can try turning off the restraints or use
only positionally independent restraints (like h-bonds) on the initial
structures. Or only restrain one of the molecules (i.e. protein).

Now, looking at the RMSd plots, we must wonder if differences of 1-3
Angstrom are significant? Atoms have net thermal motion so we expect some
RMSd; also, the more atoms, the larger the RMSd will be. To check this
out you need to look at average structures and compare them from the
different simulations. I suggest creating average structures from 3-4ns
and comparing them. You may also consider looking at 2D-RMS plots (either
across an individual trajectory or all pooled together to show
similarities).

Finally, as pointed out by Dan Roe, 4 ns is clearly not sufficient to
relax the protein, DNA, solvent and ion environment. This may require
longer...

Note also:

> Heating Stage 1
> &cntrl
> cut = 12.0,
> /

The cutoff is rather large for an explicit water simulation with PME; I
would suggest 8.0 or 9.0 angstroms, but first convince yourself that the
small differences in the results between modest changes in time step or
Langevin coupling are real.

[Instabilities could also result from the force field; I would suggest
ff99SB + parmbsc0 for the protein-DNA].

--tec3


_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Jul 08 2011 - 15:00:04 PDT
Custom Search