[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat
Bernhard
b.reuter at uni-kassel.de
Wed Jul 15 17:57:35 CEST 2015
I now also tried with bigger manual rlist values of 1.02 and 1.05: This
doesnt help - the drift of the conserved energy with Nose-Hoover gets
even bigger: -1,8 kJ/mol/ps (7.9*10^-5 kJ/mol/ps per atom) and -2.62
kJ/mol/ps (1.15*10^-4 kJ/mol/ps per atom).
Am 15/07/15 um 17:41 schrieb Bernhard:
> 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.
>>
>
More information about the gromacs.org_gmx-developers
mailing list