[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