[gmx-users] nve energy conservation

Mon Mar 19 22:25:15 CET 2018

I just want to clarify my questions:

1. Does running the simulation with NO 'settles' and NO 'constraints'
defined in the topology file, just 'bonds' and 'angles', give an accurate
system of rigid SPCE water?

2. I am not sure how to define 'bonds' and 'angles' in the topology file
that is rigid - I assume the force constant should be 0 for rigid
molecules.  However, when I set the force constant to 0 for bonds and
angles, the simulation gives a segmentation fault.

3. It appears that having constraints explicitly defined in the topology
file results in no energy conservation. Can someone explain to me why this
might be?

Below is my topology file:

[ defaults ]
; nbfunc    comb-rule    gen-pairs    fudgeLJ    fudgeQQ
  1             2             no         1.0     1.0

[ atomtypes ]
OW            15.99940      -0.8476   A           0.3166          0.65
HW            1.00800       0.4238    A           0.0             0.0

[moleculetype]
; name nrexcl
SOL 1

[atoms]
; nr type resnr residu atom cgnr charge mass
     1  OW   1    SOL    OW1      1      -0.8476     15.99940
     2  HW   1    SOL    HW2      1      0.4238     1.00800
     3  HW   1    SOL    HW3      1      0.4238     1.00800

;[ settles ]
; OW    funct   doh     dhh
;1       1       0.1     0.163298

; [constraints]
; i     j       funct   length
;1     2       2       0.1
;1     3       2       0.1

 [ bonds ]
; i     j       funct   length  force.c.
1       2       1       0.1      345000 0.1     345000
1       3       1       0.1      345000 0.1     345000

[ angles ]
; i     j       k       funct   angle   force.c.
2       1       3       1       109.47  383 109.47  383

[ exclusions ]
1       2       3
2       1       3
3       1       2

[ system ]
SPCE water

[ molecules ]
SOL               1000

On Mon, Mar 19, 2018 at 5: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.
> 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.
>>
>
>