[gmx-users] nve energy conservation

Justin Lemkul jalemkul at vt.edu
Mon Mar 19 23:57:31 CET 2018



On 3/19/18 5:25 PM, Jo wrote:
> 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?

No, exactly the opposite. Bonds and angles are modeled as harmonic 
potentials; by definition they are not rigid.

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

You don't. If you define a zero force constant for a bonded 
interactions, that means that interaction has no effect. Your molecules 
are just splitting apart if you do this, hence the seg 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?

Mark's results from a previous message dispute the assertion of "no 
energy conservation." You're always going to observe some drift, but it 
should be small. Does your GROMACS installation pass all regression 
tests? I haven't been paying much attention to this thread, but are you 
using the latest version of GROMACS with appropriate .mdp options?

-Justin

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

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.biochem.vt.edu/people/faculty/JustinLemkul.html

==================================================



More information about the gromacs.org_gmx-users mailing list