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

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