[gmx-users] nve energy conservation

[Jo]

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.
Previously, the total energy of the system was consistently increasing or
decreasing by a significant amount resulting in no energy conservation.

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.
However, I don't understand why or how removing the constraints allowed for
energy conservation. Does simply defining 'bonds' and 'angles' actually
force the atoms to to constrained as a rigid molecule.  If not, why does
adding some sort of constraint make the system loose energy conservation?
How can I go about constraining the atoms as a rigid water molecule so that
I don't loose energy conservation?

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