[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat
Bernhard
b.reuter at uni-kassel.de
Wed Jul 15 18:06:44 CEST 2015
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