[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat
Bernhard
b.reuter at uni-kassel.de
Thu Jul 16 11:37:43 CEST 2015
I see... I didn't consider the implicit constraints in rigid water
models like TIP3P (which I used). My bad! Thx for pointing that out, Mark!
I rerun the simulations (NVT, Nose-Hoover) in double precision and still
get perfect linear drifts of the Conserved Energy - but upwards and with
different reaction to variations of rlist.
To compare (rlist= auto means verlet-buffer-drift=0.005; all energies in
kJ/mol/ps per atom):
rlist Conserved-En (SP) Conserved-En (DP)
auto -7.14*10^-5 4.31*10^-5
1.012 -2.94*10^-5 8.82*10^-5
1.02 -7.94*10^-5 4.27*10^-5
1.05 -1.15*10^-4 4.93*10^-7
Interesting is the reversion of the slope for DP compared to SP.
Even more interesting is the two orders of magnitude smaller drift of
4.93*10^-7 for rlist=1.05 for DP compared to the average (even three
orders smaller compared to the SP value of -1.15*10^-4 for rlist=1.05).
Since the drift of the conserved energy for DP and rlist=1.05 is only
5*10^-4 % of the average conserved energy in 100ps simulation time, one
could even say that the conserved energy is constant in this case.
While in the SP case the drift for rlist=1.05 is the biggest and quite
significant.
Seems like the nearly perfect conservation for DP and rlist=1.05 is
something like a lucky cancellation of errors?
My conclusions of this are at least:
- The drift isn't resulting from a round-off error (since DP doesn't
improve the behaviour significantly).
- It may be an algorithmic error (since it goes linear with the number
of steps N) that is at least influenced by the size of the verlet
buffer (rlist).
- Its not related to the thermostat since NVE doesn't change anything.
But were is the reversion of the slope for DP compared to SP coming from?
Can this really be related to the constraints (mainly of the water since
switching of the constraints for the protein doesn't change anything)?
At last a simple technical question: Were can I find the value of rlist
if verlet-buffer-drift is switched on and rlist is calculated automatically?
Best,
Bernhard
Am 16/07/15 um 08:10 schrieb Mark Abraham:
>
> Hi,
>
> Likely your water still used SETTLE, as the mdp option merely disables
> the automatic conversion of normal bonds to constraints.
>
> Mark
>
>
> On Wed, 15 Jul 2015 19:12 Bernhard <b.reuter at uni-kassel.de
> <mailto:b.reuter at uni-kassel.de>> wrote:
>
> Finally I also tried a NVE simulation (rlist=1.012nm) with the
> same system, parameters and starting conditions but _without__any
> constraints_ (constraints=none): The result is an equal linear
> drift (-0.78 kJ/mol/ps or -3.43*10^-5 kJ/mol/ps per atom) of the
> total energy as for NVE with constraints=h-bonds or of the
> conserved energy for Nose-Hoover with constraints=h-bonds
> (rlist=1.012nm in all cases; lincs-iter=2 for NVE and lincs-iter=1
> for NVT).
> So constraints should be ruled out as cause of the drift...
>
> Best,
> Bernhard
>
> Am 15/07/15 um 18:14 schrieb Mark Abraham:
>> 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>.
>>
>>
>>
>
> --
> 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/20150716/6bb428d2/attachment-0001.html>
More information about the gromacs.org_gmx-developers
mailing list