[gmx-users] nve energy conservation

[Mark Abraham]

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.
>