[gmx-users] nve energy conservation

Tue Mar 20 21:02:41 CET 2018

Thanks everyone for your reply.
I see your point that it could just be some error cancellation from
constraints and parameters.

I was previously using gromacs-5.1.4, and changed to gromacs-2018, and the
drift seems to be much smaller! I get something similar to what Mark had
posted -0.000025 kJ/mol/ps/atom for 1000 SPCE waters with settles at 298K
for single precision.  With double precision, there is no noticeable drift
for a short test run of only 200 ps. Any idea why this is? Is there some
default setting that changed?

Mark, regarding your earlier questions: My application is to calculate very
precise and accurate dynamic properties (diffusion, viscosity) of NaCl in
water. There are a number of works in literature where people use NVT
simulations to calculate these properties, but I wanted to check if the
thermostat affected the dynamics any by running an NVE simulation.  Also I
want to better understand how NVE simulations work and all the controlling
factors.

>From figure 3.5 in the reference manual, I think 10^-5 drift (kJ/mol/ps) is
acceptable.

Can you help me better understand what verlet-buffer-tolerance does?  My
understanding from the documentation is that it determines the width of the
buffer outside of the cutoff for the neighborlist, but I am not clear why
energy can be added or lost by manipulating this parameter.

Also, I am not sure how to judge what is an acceptable width of fluctuation
NVE simulations?  Does it correspond to random error (from the paper
<https://www.biorxiv.org/content/biorxiv/early/2016/10/24/083055.full.pdf>Mark
referenced)?

Thank you!

On Tue, Mar 20, 2018 at 1:39 AM, Mark Abraham <mark.j.abraham at gmail.com>
wrote:

> Hi,
>
> On Mon, Mar 19, 2018 at 10:08 PM Jo <jojo412202 at gmail.com> wrote:
>
> > Hello,
> >
> > Thank you for your response.  I have found a way to conserve energy of my
> > box of SPC/E water by removing 'settles' from my topology file.
> >
>
> Thus your water is no longer rigid.
>
>
> > Previously, the total energy of the system was consistently increasing or
> > decreasing by a significant amount resulting in no energy conservation.
> >
>
> You're doing a numerical computation with finite-precision forces,
> positions and velocities and a discrete fixed time step. You are guaranteed
> to have both random and systematic errors compared to exact arithmetic.
> Maybe you can find a parameter combination that will appear to lack
> systematic errors, but that doesn't mean errors are absent, or the
> combination is "better," as you have no doubt already read about in the
> article I shared.
>
>
> > I have defined 'bonds' and 'angles' with length and force constant for a
> > flexible SPC/E water, with no 'settles' and no 'constraints'.  For these
> > parameters, the total energy converges to approximately the correct
> energy
> > and fluctuates withing 40 kJ/mol for a timestep of 1 fs, which is good.
> >
>
> Why is it good? :-)
>
>
> > However, I don't understand why or how removing the constraints allowed
> for
> > energy conservation.
>
>
> It is implementing different equations of motion. If the only observable of
> interest is the potential energy, and the multiple sources of numerical
> inaccuracy are cancelling suitable for you, and you can afford to generate
> enough samples at that short time step, then go right ahead :-)
>
>
> > Does simply defining 'bonds' and 'angles' actually
> > force the atoms to to constrained as a rigid molecule.
>
>
> No, see the documentation about bonded interactions...
>
>
> >   If not, why does
> > adding some sort of constraint make the system loose energy conservation?
> >
>
> It is implementing different equations of motion, and they have different
> numerical properties. Whether those properties are sufficient for any
> purpose is an open question. I note that you ignored my earlier series of
> questions on this topic, which were intended to help you learn ;-)
>
> How can I go about constraining the atoms as a rigid water molecule so that
> > I don't loose energy conservation?
> >
>
> As you can infer from Figure 3.5 in the reference manual, you can choose
> the drift from the verlet buffer to add enough energy to offset the drift
> from SETTLE. But obviously that would be a pair of compensating systematic
> errors - so you won't learn anything from merely observing
> that the systematic drift is zero.
>
> Mark
>
>
> > Thank you in advance,
> >
> > Jo
> >
> > On Fri, Mar 16, 2018 at 4:51 PM, Mark Abraham <mark.j.abraham at gmail.com>
> > wrote:
> >
> > > Hi,
> > >
> > > On Fri, Mar 16, 2018 at 7:42 PM Jo <jojo412202 at gmail.com> wrote:
> > >
> > > > Thank you for your reply!
> > > >
> > > > I am attempting to conserve energy in an NVE run of 1000 SPCE water
> - I
> > > > have tried a number of different verlet-buffer-tolerances (0.001 to
> > > 5e-5),
> > > > sometimes the run output file suggests a specific verlet
> > > buffer-tolerance.
> > > > However I am still experiencing ~600 kJ/mol shift per ns.
> > >
> > >
> > > You can't be getting that over that whole range. If you're doing NVE
> with
> > > tolerance 5e-5 then the drift is negative (from SETTLE, e.g. I just
> > > observed -0.000069 kJ/mol/ps/atom on 1728 tip3p waters, GROMACS 2018
> > mixed
> > > precision on a GPU, with SETTLE at 300K and 2fs timestep with PME and
> > other
> > > settings at defaults).
> > >
> > > I have tried
> > > > double precision which does not seem to make a difference.  I also
> use
> > a
> > > > potential modifier (potential-shift) to ensure the potential reaches
> 0.
> > >
> > >
> > > You can never have an unshifted potential with the Verlet scheme.
> > >
> > >
> > > > I
> > > > have tried turning off the charges (to see if the ewald parameters
> are
> > > the
> > > > source of the problem), but the same massive energy shift occurs -
> > which
> > > > leads me to believe the ewald parameters are not the source of the
> > > problem.
> > > >
> > > > I am using 'settle' to constrain the water molecule.  I am suspicious
> > > that
> > > > this could be the cause of the energy shift. Although, looking at
> some
> > of
> > > > the previous posts on gromacs email forums, it appears that 'settle'
> > has
> > > > not been a problem for energy conservation previously.
> > > >
> > > > I have also tried different versions of Gromacs (5.1.4 and 2018), but
> > the
> > > > energy shift still occurs.
> > > >
> > > > Can you recommend any other parameters to change?  Or any other
> > > strategies
> > > > to go about conserving energy for NVE?
> > >
> > >
> > > First, given that you can't have zero drift in a numerical simulation
> > with
> > > a finite time step (and particularly not with constraints), have you
> > > decided what level you regard as acceptable? For what application? Have
> > you
> > > considered the thoughts at https://www.biorxiv.org/node/23099.full?
> > >
> > > Mark
> > > --
> > > Gromacs Users mailing list
> > >
> > > * Please search the archive at http://www.gromacs.org/
> > > Support/Mailing_Lists/GMX-Users_List before posting!
> > >
> > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> > >
> > > * For (un)subscribe requests visit
> > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > > send a mail to gmx-users-request at gromacs.org.
> > >
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > send a mail to gmx-users-request at gromacs.org.
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>