[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat
Bernhard
b.reuter at uni-kassel.de
Wed Jul 15 18:16:26 CEST 2015
I remember some research article http://dx.doi.org/10.1063/1.2431176 (or
https://www.deshawresearch.com/publications/A%20common,%20avoidable%20source%20of%20error%20in%20molecular%20dynamics%20integrators.pdf)
from 2006 which showed a comparable linear energy drift for single
precision GROMACS 3.3.1 due to not optimal calculation of the velocity
of constrained particels.
But this issue should be solved in version 4.6.7?
Best,
Bernhard
Am 15/07/15 um 18:04 schrieb Bernhard:
> Indeed seems so... unfortunately I have no clue about the cause.
> Maybe the head of the .log file is of some use?
>
> Log file opened on Wed Jul 1 17:33:48 2015
> Host: theo2-pc20 pid: 15411 nodeid: 0 nnodes: 1
> Gromacs version: VERSION 4.6.7
> Precision: single
> Memory model: 64 bit
> MPI library: thread_mpi
> OpenMP support: enabled
> GPU support: disabled
> invsqrt routine: gmx_software_invsqrt(x)
> CPU acceleration: AVX_256
> FFT library: fftw-3.3.2-sse2
> Large file support: enabled
> RDTSCP usage: enabled
> Built on: Di 16. Jun 15:38:40 CEST 2015
> Built by: berni at theo2-pc20 [CMAKE]
> Build OS/arch: Linux 3.13.0-43-generic x86_64
> Build CPU vendor: GenuineIntel
> Build CPU brand: Intel(R) Core(TM) i7-4930K CPU @ 3.40GHz
> Build CPU family: 6 Model: 62 Stepping: 4
> Build CPU features: aes apic avx clfsh cmov cx8 cx16 f16c htt lahf_lm
> mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdrnd rdtscp
> sse2 sse3 sse4.1 sse4.2 ssse3 tdt x2apic
> C compiler: /usr/bin/cc GNU cc (Ubuntu 4.8.2-19ubuntu1) 4.8.2
> C compiler flags: -mavx -Wextra -Wno-missing-field-initializers
> -Wno-sign-compare -Wall -Wno-unused -Wunused-value
> -Wno-unused-parameter -Wno-array-bounds -Wno-maybe-uninitialized
> -Wno-strict-overflow -fomit-frame-pointer -funroll-all-loops
> -fexcess-precision=fast -O3 -DNDEBUG
>
>
> :-) G R O M A C S (-:
>
> GROwing Monsters And Cloning Shrimps
>
> :-) VERSION 4.6.7 (-:
>
> Contributions from Mark Abraham, Emile Apol, Rossen Apostolov,
> Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,
> Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph
> Junghans,
> Peter Kasson, Carsten Kutzner, Per Larsson, Pieter Meulenhoff,
> Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz,
> Michael Shirts, Alfons Sijbers, Peter Tieleman,
>
> Berk Hess, David van der Spoel, and Erik Lindahl.
>
> Copyright (c) 1991-2000, University of Groningen, The Netherlands.
> Copyright (c) 2001-2012,2013, The GROMACS development team at
> Uppsala University & The Royal Institute of Technology, Sweden.
> check out http://www.gromacs.org for more information.
>
> This program is free software; you can redistribute it and/or
> modify it under the terms of the GNU Lesser General Public License
> as published by the Free Software Foundation; either version 2.1
> of the License, or (at your option) any later version.
>
> :-) mdrun (-:
>
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> B. Hess and C. Kutzner and D. van der Spoel and E. Lindahl
> GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable
> molecular simulation
> J. Chem. Theory Comput. 4 (2008) pp. 435-447
> -------- -------- --- Thank You --- -------- --------
>
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H.
> J. C.
> Berendsen
> GROMACS: Fast, Flexible and Free
> J. Comp. Chem. 26 (2005) pp. 1701-1719
> -------- -------- --- Thank You --- -------- --------
>
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> E. Lindahl and B. Hess and D. van der Spoel
> GROMACS 3.0: A package for molecular simulation and trajectory analysis
> J. Mol. Mod. 7 (2001) pp. 306-317
> -------- -------- --- Thank You --- -------- --------
>
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> H. J. C. Berendsen, D. van der Spoel and R. van Drunen
> GROMACS: A message-passing parallel molecular dynamics implementation
> Comp. Phys. Comm. 91 (1995) pp. 43-56
> -------- -------- --- Thank You --- -------- --------
>
> Input Parameters:
> integrator = md
> nsteps = 10000000
> init-step = 0
> cutoff-scheme = Verlet
> ns_type = Grid
> nstlist = 10
> ndelta = 2
> nstcomm = 100
> comm-mode = Linear
> nstlog = 5000
> nstxout = 5000
> nstvout = 5000
> nstfout = 0
> nstcalcenergy = 100
> nstenergy = 1000
> nstxtcout = 1000
> init-t = 0
> delta-t = 0.001
> xtcprec = 1000
> fourierspacing = 0.12
> nkx = 60
> nky = 60
> nkz = 60
> pme-order = 4
> ewald-rtol = 1e-05
> ewald-geometry = 0
> epsilon-surface = 0
> optimize-fft = FALSE
> ePBC = xyz
> bPeriodicMols = FALSE
> bContinuation = TRUE
> bShakeSOR = FALSE
> etc = Nose-Hoover
> bPrintNHChains = FALSE
> nsttcouple = 10
> epc = No
> epctype = Isotropic
> nstpcouple = -1
> tau-p = 1
> ref-p (3x3):
> ref-p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> ref-p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> compress (3x3):
> compress[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> compress[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> compress[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> refcoord-scaling = No
> posres-com (3):
> posres-com[0]= 0.00000e+00
> posres-com[1]= 0.00000e+00
> posres-com[2]= 0.00000e+00
> posres-comB (3):
> posres-comB[0]= 0.00000e+00
> posres-comB[1]= 0.00000e+00
> posres-comB[2]= 0.00000e+00
> verlet-buffer-drift = 0.005
> rlist = 1
> rlistlong = 1
> nstcalclr = 10
> rtpi = 0.05
> coulombtype = PME
> coulomb-modifier = Potential-shift
> rcoulomb-switch = 0
> rcoulomb = 1
> vdwtype = Cut-off
> vdw-modifier = Potential-shift
> rvdw-switch = 0
> rvdw = 1
> epsilon-r = 1
> epsilon-rf = inf
> tabext = 1
> implicit-solvent = No
> gb-algorithm = Still
> gb-epsilon-solvent = 80
> nstgbradii = 1
> rgbradii = 1
> gb-saltconc = 0
> gb-obc-alpha = 1
> gb-obc-beta = 0.8
> gb-obc-gamma = 4.85
> gb-dielectric-offset = 0.009
> sa-algorithm = Ace-approximation
> sa-surface-tension = 2.05016
> DispCorr = EnerPres
> bSimTemp = FALSE
> free-energy = no
> nwall = 0
> wall-type = 9-3
> wall-atomtype[0] = -1
> wall-atomtype[1] = -1
> wall-density[0] = 0
> wall-density[1] = 0
> wall-ewald-zfac = 3
> pull = no
> rotation = FALSE
> disre = No
> disre-weighting = Conservative
> disre-mixed = FALSE
> dr-fc = 1000
> dr-tau = 0
> nstdisreout = 100
> orires-fc = 0
> orires-tau = 0
> nstorireout = 100
> dihre-fc = 0
> em-stepsize = 0.01
> em-tol = 10
> niter = 20
> fc-stepsize = 0
> nstcgsteep = 1000
> nbfgscorr = 10
> ConstAlg = Lincs
> shake-tol = 0.0001
> lincs-order = 4
> lincs-warnangle = 30
> lincs-iter = 1
> bd-fric = 0
> ld-seed = 1993
> cos-accel = 0
> deform (3x3):
> deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
> adress = FALSE
> userint1 = 0
> userint2 = 0
> userint3 = 0
> userint4 = 0
> userreal1 = 0
> userreal2 = 0
> userreal3 = 0
> userreal4 = 0
> grpopts:
> nrdf: 1501.9 44334.1
> ref-t: 300 300
> tau-t: 2.5 2.5
> anneal: No No
> ann-npoints: 0 0
> acc: 0 0 0
> nfreeze: N N N
> energygrp-flags[ 0]: 0
> efield-x:
> n = 0
> efield-xt:
> n = 0
> efield-y:
> n = 0
> efield-yt:
> n = 0
> efield-z:
> n = 0
> efield-zt:
> n = 0
> bQMMM = FALSE
> QMconstraints = 0
> QMMMscheme = 0
> scalefactor = 1
> qm-opts:
> ngQM = 0
> Using 1 MPI thread
> Using 12 OpenMP threads
>
> Detecting CPU-specific acceleration.
> Present hardware specification:
> Vendor: GenuineIntel
> Brand: Intel(R) Core(TM) i7-4930K CPU @ 3.40GHz
> Family: 6 Model: 62 Stepping: 4
> Features: aes apic avx clfsh cmov cx8 cx16 f16c htt lahf_lm mmx msr
> nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdrnd rdtscp sse2
> sse3 sse4.1 sse4.2 ssse3 tdt x2apic
> Acceleration most likely to fit this hardware: AVX_256
> Acceleration selected at GROMACS compile time: AVX_256
>
> Will do PME sum in reciprocal space.
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G.
> Pedersen
> A smooth particle mesh Ewald method
> J. Chem. Phys. 103 (1995) pp. 8577-8592
> -------- -------- --- Thank You --- -------- --------
>
> Will do ordinary reciprocal space Ewald sum.
> Using a Gaussian width (1/beta) of 0.320163 nm for Ewald
> Cut-off's: NS: 1 Coulomb: 1 LJ: 1
> Long Range LJ corr.: <C6> 2.8855e-04
> System total charge: 0.000
> Generated table with 1000 data points for Ewald.
> Tabscale = 500 points/nm
> Generated table with 1000 data points for LJ6.
> Tabscale = 500 points/nm
> Generated table with 1000 data points for LJ12.
> Tabscale = 500 points/nm
> Generated table with 1000 data points for 1-4 COUL.
> Tabscale = 500 points/nm
> Generated table with 1000 data points for 1-4 LJ6.
> Tabscale = 500 points/nm
> Generated table with 1000 data points for 1-4 LJ12.
> Tabscale = 500 points/nm
>
> Using AVX-256 4x4 non-bonded kernels
>
> Using Lorentz-Berthelot Lennard-Jones combination rule
>
> Potential shift: LJ r^-12: 1.000 r^-6 1.000, Ewald 1.000e-05
> Initialized non-bonded Ewald correction tables, spacing: 6.60e-04
> size: 3033
>
> Pinning threads with an auto-selected logical core stride of 1
>
> Initializing LINear Constraint Solver
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> B. Hess and H. Bekker and H. J. C. Berendsen and J. G. E. M. Fraaije
> LINCS: A Linear Constraint Solver for molecular simulations
> J. Comp. Chem. 18 (1997) pp. 1463-1472
> -------- -------- --- Thank You --- -------- --------
>
> The number of constraints is 292
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> S. Miyamoto and P. A. Kollman
> SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithms for
> Rigid
> Water Models
> J. Comp. Chem. 13 (1992) pp. 952-962
> -------- -------- --- Thank You --- -------- --------
>
> Center of mass motion removal mode is Linear
> We have the following groups for center of mass motion removal:
> 0: rest
> There are: 22765 Atoms
> Initial temperature: 299.915 K
>
>
> Am 15/07/15 um 17:57 schrieb Shirts, Michael R. (mrs5pt):
>>> 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...
>>
>> So sounds like the conserved quantity drift with NVT is not due to
>> NH, but
>> due to something else with the underlying dynamics.
>>
>> Best,
>> ~~~~~~~~~~~~
>> Michael Shirts
>> Associate Professor
>> Department of Chemical Engineering
>> University of Virginia
>> michael.shirts at virginia.edu
>> (434) 243-1821
>>
>>
>>
>> On 7/15/15, 11:41 AM, "Bernhard" <b.reuter at uni-kassel.de> wrote:
>>
>>> 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.
>>> --
>>> 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