On Tue, Apr 7, 2015 at 3:52 PM, Steven Ramsey <vpsramsey.gmail.com> wrote:
> Dear Amber Users,
>
> Sorry for the potential resubmission,
>
> I am attempting to simulate a protein in NVT, while saving the
> configuration every 5ps using the langevin thermostat. My problem is I
> found a drift in average total energy and am uncertain as to whether this
> drift is due to some error on my part or solely due to the behavior of the
> thermostat.
>
> These ETOT values were collected by process_mdout.perl from the .out files
> of 5 contiguous simulations of 20k frames each run via the input file
> attached (md10.in).
>
>
> Also attached is the xmgrace generated plot of the running ETOT average.
>
>
>
> Does the total energy reported in the .out files include the restraint
>
> energies (my system is heavily restrained on protein heavy atoms) and would
>
> this contribute to this drift? Is there any blaring errors in my .in file
>
> or my protocol (100ns while saving configurations every 5ps done 5x) which
>
> could result in this drift?
>
It certainly can. Energy drift in NVT simulations often indicates that a
simulation is evolving. (Note that random fluctuations are normal for a
canonical ensemble, but a systematic drift is not). I would interpret a
systematic drift as the system equilibrating (because if it's sampling from
the surface of its phase space, then it would have only random fluctuations
about its mean).
Some things to check: if you are running NVT and the thermostat is doing
its job, the kinetic energy should be fairly constant (*that* should
certainly not be evolving). Which means that the potential energy -- which
is *strictly* a function of the protein's conformation -- is contributing
exclusively to the drift. Which means that conformational changes are the
root cause of this drift.
Note that with restraints, the free energy surface is different from that
of the unrestrained system (since the restraints are now a part of the
potential energy function). So it's definitely possible for the restraints
to be the culprit here -- plot the restraint energy to see if there is a
trend to that.
Also visualize your system -- see if you can detect this irreversible
change that your system is undergoing just by looking. That could help you
determine what's happening.
Also, I would suggest having a look at mdout_analyzer.py instead of
process_mdout.pl (just install matplotlib and numpy first). It's a
simpler, more interactive GUI tool for investigating energy data from MD
simulations.
> For references sake the system was equilibrated in NPT conditions prior to
>
> these production NVT simulations for ~1.4ns and the forcefields in use are
>
> TIP3P for the water and AMBER99SB for the protein.
>
This is a short equilibration. Based on your results, I would say it's
clearly not long enough.
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Apr 07 2015 - 13:30:03 PDT
