[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat
Bernhard
b.reuter at uni-kassel.de
Wed Jul 15 21:18:03 CEST 2015
VDW is shifted (vdw-modifier=Potential-shift). The simulations are run
in single precision.
Currently Im re-runing them in DP.
Best,
Bernhard
Am 15/07/15 um 21:09 schrieb Shirts, Michael R. (mrs5pt):
> Just wondering, Is this with single or double precision?
>
> Also, is the VDW interaction tapered or shifted? (check the log file
> or mdp out) If not, that can easily cause lack of energy conservation.
>
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Associate Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu
> (434) 243-1821
>
> From: Bernhard <b.reuter at uni-kassel.de <mailto:b.reuter at uni-kassel.de>>
> Reply-To: "gmx-developers at gromacs.org
> <mailto:gmx-developers at gromacs.org>" <gmx-developers at gromacs.org
> <mailto:gmx-developers at gromacs.org>>
> Date: Wednesday, July 15, 2015 at 11:09 AM
> To: "gromacs.org_gmx-developers at maillist.sys.kth.se
> <mailto:gromacs.org_gmx-developers at maillist.sys.kth.se>"
> <gromacs.org_gmx-developers at maillist.sys.kth.se
> <mailto:gromacs.org_gmx-developers at maillist.sys.kth.se>>
> Subject: Re: [gmx-developers] Drift in Conserved-Energy with
> Nose-Hoover thermostat
>
> 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
>> <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
>> <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
>> <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/20150715/0ce0214f/attachment-0001.html>
More information about the gromacs.org_gmx-developers
mailing list