[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat

[Shirts, Michael R. (mrs5pt)]

> There I also got a linear drift (but smaller) of 0.78 kJ/mol/ps
>(3.436*10^-5 kJ/mol/ps per atom).

> For comparison reasons I also did a NVT Nose-Hoover Simulation with
>manually set rlist=1.012nm: There I got a comparable linear drift of 0.67
>kJ/mol/ps (2.94*10^-5 kJ/mol/ps per atom). So no differences between NVE
>and NVT so far in my opinion...


So sounds like the conserved quantity drift with NVT is not due to NH, but
due to something else with the underlying dynamics.

Best,
~~~~~~~~~~~~
Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434) 243-1821



On 7/15/15, 11:41 AM, "Bernhard" <b.reuter at uni-kassel.de> wrote:

>I mean a drift of the total energy in NVE - while with Nose-Hoover the
>drift is in the Conserved-Energy quantity of g_energy (the total energy
>shows no drift with Noose-Hoover...).
>
>Am 15/07/15 um 17:38 schrieb Bernhard:
>> I also did a NVE simulation with the same parameters, system and
>> starting conditions but with manually set rlist=1.012nm (since
>> verlet-buffer-drift doesnt work in NVE):
>>There I also got a linear drift (but smaller) of 0.78 kJ/mol/ps
>>(3.436*10^-5 kJ/mol/ps per atom).
>>For comparison reasons I also did a NVT Nose-Hoover Simulation with
>>manually set rlist=1.012nm:
>>There I got a comparable linear drift of 0.67 kJ/mol/ps (2.94*10^-5
>>kJ/mol/ps per atom).
>>So no differences between NVE and NVT so far in my opinion...
>>
>>
>>
>> Best,
>> Bernhard
>>
>> Am 15/07/15 um 17:10 schrieb Shirts, Michael R. (mrs5pt):
>>> The conserved quantity in nose-hoover is not quite as good as the
>>> conserved energy, which should have no drift at all.  For NH, the
>>> conserved quantity should drift as a random Gaussian process with mean
>>> zero (i.e. go with sqrt(N)). It shouldn't be drifting linearly.
>>>
>>> I would check to see if your system conserved energy when run with NVE
>>> (use the endpoint of the NPT simulation).  It's easier to diagnose any
>>> problems with an NVE simulation, which should have virtually no
>>> drift, vs
>>> a NVT simulation, which has random noise drift.  Odds are, if there is
>>>a
>>> problem with the NVT simulation, it will also show up in the NVE
>>> simulation if only the thermostat is removed.
>>>
>>> Also, consider looking at http://pubs.acs.org/doi/abs/10.1021/ct300688p
>>> for tests of whether the ensemble generated is correct.
>>>
>>> Best,
>>> ~~~~~~~~~~~~
>>> Michael Shirts
>>> Associate Professor
>>> Department of Chemical Engineering
>>> University of Virginia
>>> michael.shirts at virginia.edu
>>> (434) 243-1821
>>>
>>>
>>>
>>> On 7/15/15, 10:58 AM, "Bernhard" <b.reuter at uni-kassel.de> wrote:
>>>
>>>> Dear Gromacs Users and Developers,
>>>>
>>>> I have a problem regarding energy conservation in my 10ns NVT
>>>> protein+water+ions (22765 atoms) production (minimization and
>>>> equilibration for more than 15ns was carried out in NPT before)
>>>> simulations using a Nose-Hoover thermostat (tau=2.5ps).
>>>> On first glance everything looks fine - the potential, kinetic and
>>>> total
>>>> energy are nearly perfectly constant (with normal fluctuations) - but
>>>> when I checked the "Conserved-Energy" quantity that g_energy outputs I
>>>> had to recognize a significant (nearly perfectly) linear downward
>>>>drift
>>>> of this "to-be-conserved" quantity of around 1.7 kJ/mol/ps (7.48*10^-5
>>>> kJ/mol/ps per atom).
>>>> This appears somehow disturbing to me since I would expect that this
>>>> Conserved-Energy is the conserved energy of the extended Nose-Hoover
>>>> Hamiltonian - which should by definition be conserved.
>>>>
>>>> If it would be a drift caused by normal round-off error due to single
>>>> precision I would expect it to grow with Sqrt(N) and not with N
>>>> (linear)
>>>> (N=number of steps).
>>>> So I would like to know if this is a normal behaviour and also what
>>>> could cause this (buffer size, precision, constraints etc)?
>>>> Also I would like to know, if I am correct with my guess that the
>>>> "Conserved-Energy" quantity is in this case the energy of the extended
>>>> Nose-Hoover Hamiltonian?
>>>> The .mdp file is atatched (don't be confused about rlist=1 - since Im
>>>> using the Verlet-scheme the verlet-buffer-drift option should be by
>>>> default active and determine the rlist value (Verlet buffer-size)
>>>> automatically).
>>>>
>>>> Best regards,
>>>> Bernhard
>>>>
>>>> ; Run parameters
>>>> integrator    = md        ; leap-frog integrator
>>>> nsteps        = 10000000    ; 10000 ps = 10 ns
>>>> dt        = 0.001        ; 1 fs
>>>> ; Output control
>>>> nstxout        = 5000        ; save coordinates every ps
>>>> nstvout        = 5000        ; save velocities every ps
>>>> nstxtcout    = 1000        ; xtc compressed trajectory output every ps
>>>> nstenergy    = 1000        ; save energies every ps
>>>> nstlog        = 5000        ; update log file every ps
>>>> ; Bond parameters
>>>> continuation    = yes        ; continue from NPT
>>>> constraint_algorithm = lincs    ; holonomic constraints
>>>> constraints    = h-bonds    ; all bonds (even heavy atom-H bonds)
>>>> constrained
>>>> lincs_iter    = 1        ; accuracy of LINCS
>>>> lincs_order    = 4        ; also related to accuracy
>>>> ; Neighborsearching
>>>> cutoff-scheme    = Verlet    ; Verlet cutoff-scheme instead of
>>>> group-scheme (no charge-groups used)
>>>> ns_type        = grid        ; search neighboring grid cells
>>>> nstlist        = 10        ; 10 fs
>>>> rlist        = 1.0        ; short-range neighborlist cutoff (in nm)
>>>> rcoulomb    = 1.0        ; short-range electrostatic cutoff (in nm)
>>>> rvdw        = 1.0        ; short-range van der Waals cutoff (in nm)
>>>> ; Electrostatics
>>>> coulombtype    = PME        ; Particle Mesh Ewald for long-range
>>>> electrostatics
>>>> pme_order    = 4        ; cubic interpolation
>>>> fourierspacing    = 0.12        ; grid spacing for FFT
>>>> ; Temperature coupling is on
>>>> tcoupl        = nose-hoover    ; modified Berendsen thermostat
>>>> tc-grps        = Protein Non-Protein    ; two coupling groups - more
>>>> accurate
>>>> tau_t        = 2.5    2.5    ; time constant, in ps
>>>> ref_t        = 300     300    ; reference temperature, one for each
>>>> group, in K
>>>> ; Pressure coupling is off
>>>> pcoupl        = no         ; no pressure coupling in NVT
>>>> ; Periodic boundary conditions
>>>> pbc        = xyz        ; 3-D PBC
>>>> ; Dispersion correction
>>>> DispCorr    = EnerPres    ; account for cut-off vdW scheme
>>>> ; Velocity generation
>>>> gen_vel        = no        ; don¹t assign velocities from Maxwell
>>>> distribution
>>>>
>>>> -- 
>>>> Gromacs Developers mailing list
>>>>
>>>> * Please search the archive at
>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers
>>>> or send a mail to gmx-developers-request at gromacs.org.
>>
>
>-- 
>Gromacs Developers mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers
>or send a mail to gmx-developers-request at gromacs.org.