[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat
Bernhard Reuter
b.reuter at uni-kassel.de
Sun Jul 19 22:13:30 CEST 2015
Dear All,
after some more literature recherche and considering the article Mark
mentioned I came to the conclusion that the systematic linear energy
drift (-0.074 kJ/mol/ns/atom or -0.03 k_B*T/ns/atom) in my NVE and NVT
simulations matches the published values: It could be considered "normal".
But - and this is a big "but" - in a 40ns simulation of my 22765 atom
protein-water-system, this drift accumulates to a 25% loss of total
energy - which is by no means negligible.
This I consider rather serious - especially considering that the md
leapfrog integrator I used for the simulations is considerd to be
symplectic, which means that it should conserve the volume in phase space.
But if you have a constant downward drift of energy you must consider
that there is less phase space volume at lower energies - so there is no
volume conservation in phase space.
This makes me doubt strongly in the physical correctness of the systems
time evolution.
So what should be the consequence if such a systematic linear energy
drift is to be considered normal for single precision GROMACS?
Since I found that in double precision the drift vanishes for a Verlet
buffer >0.05nm (while in SP it never vanishes): Do I have to switch to
double precision if I care about energy conservation, integrator
symplecticity, phase space volume conservation and ergodicity?
Since bigger round-off errors by reduced precision shouldn't accumulate
linearly but at worst with Sqrt(N): Shouldn't one be worried about the
occurence of a linear systematic error by only changing the precision
from double to single in a calculation?
Please excuse my rather critical argumentation but the mentioned points
are rather important for my research.
Please could you give me feedback on my doubts?
Best regards,
Bernhard
--
Dipl.-Phys. Bernhard Reuter
Institute of Physics
Theoretical Physics - University of Kassel
Heinrich-Plett-Str. 40
34132 Kassel - Germany
Tel.: +49-561-804-4482
Email: b.reuter at uni-kassel.de
Am 7/19/15, 9:11 PM, schrieb Bernhard Reuter:
>
>
>
> -------- Weitergeleitete Nachricht --------
> Betreff: Re: [gmx-developers] Drift in Conserved-Energy with
> Nose-Hoover thermostat
> Datum: Wed, 15 Jul 2015 16:14:40 +0000
> Von: Mark Abraham <mark.j.abraham at gmail.com>
> Antwort an: gmx-developers at gromacs.org
> An: gmx-developers at gromacs.org,
> gromacs.org_gmx-developers at maillist.sys.kth.se
>
>
>
> Hi,
>
> This looks approximately normal for constrained dynamics of systems
> with water. See e.g. fig 8 of http://dx.doi.org/10.1016/j.cpc.2013.06.003
>
> Mark
>
> On Wed, Jul 15, 2015 at 5:57 PM Bernhard <b.reuter at uni-kassel.de
> <mailto:b.reuter at uni-kassel.de>> wrote:
>
> 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 <mailto:michael.shirts at virginia.edu>
> >>> (434) 243-1821
> >>>
> >>>
> >>>
> >>> On 7/15/15, 10:58 AM, "Bernhard" <b.reuter at uni-kassel.de
> <mailto: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
> <mailto: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
> <mailto:gmx-developers-request at gromacs.org>.
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20150719/24566d77/attachment-0001.html>
More information about the gromacs.org_gmx-developers
mailing list