Project

General

Profile

Bug #1000

Possible bug with v-rescale thermostat and md-vv

Added by Mark Abraham about 5 years ago. Updated almost 4 years ago.

Status:
Closed
Priority:
Normal
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I tried the md-vv integrator in release-4-5-patches (2859895) in concert with the v-rescale (Bussi) thermostat and found some weird heating behaviour on a 900-TIP3P system that was fine with the md integrator. (.mdp file below) Is/was this thermostat known to work with md-vv?

Attached a .tpr and a tarball of input and output.

With md-vv:

Step           Time         Lambda
0 0.00000 0.00000
Energies (kJ/mol)
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
1.56788e+04 -2.52945e+02 -3.58017e+04 -3.73679e+03 -2.41127e+04
Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
9.86312e+03 -1.42495e+04 -1.42495e+04 4.39597e+02 -1.54019e+02
Pressure (bar)
-5.69569e+04
Step           Time         Lambda
500 1.00000 0.00000
Energies (kJ/mol)
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
8.23626e+03 -2.52945e+02 -1.74349e+04 -3.27515e+03 -1.27267e+04
Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
4.19512e+04 2.92245e+04 2.45247e+04 1.86976e+03 -1.54019e+02
Pressure (bar)
2.97051e+04

<snip>

Step           Time         Lambda
500000 1000.00000 0.00000

Writing checkpoint, step 500000 at Mon Sep 3 23:15:57 2012

Energies (kJ/mol)
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
1.01917e+04 -2.52945e+02 -1.51940e+04 -3.02693e+03 -8.28214e+03
Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
5.17784e+04 4.34963e+04 4.79009e+04 2.30775e+03 -1.54019e+02
Pressure (bar)
3.82536e+04
<======  ###############  >
<== A V E R A G E S ====>
<== ############### ======>
Statistics over 500001 steps using 50001 frames
Energies (kJ/mol)
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
9.71441e+03 -2.52945e+02 -1.51389e+04 -3.03171e+03 -8.70917e+03
Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
5.25444e+04 4.38352e+04 4.10799e+04 2.34189e+03 -1.54019e+02
Pressure (bar)
3.73202e+04

Obviously this is heating massively at the start and the thermostat is totally failing to deal with it. With integrator md:

Step           Time         Lambda
0        0.00000        0.00000
Energies (kJ/mol)
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      Potential
1.56788e+04   -2.52945e+02   -3.58017e+04   -3.73679e+03   -2.41127e+04
Kinetic En.   Total Energy  Conserved En.    Temperature Pres. DC (bar)
1.30540e+04   -1.10587e+04   -1.10587e+04    5.81814e+02   -1.54019e+02
Pressure (bar)
-1.29542e+04
Step           Time         Lambda
500        1.00000        0.00000
Energies (kJ/mol)
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      Potential
4.57027e+03   -2.52945e+02   -3.01897e+04   -4.38326e+03   -3.02556e+04
Kinetic En.   Total Energy  Conserved En.    Temperature Pres. DC (bar)
1.07502e+04   -1.95054e+04   -1.39359e+04    4.79135e+02   -1.54019e+02
Pressure (bar)
3.29939e+03

<snip>

Step           Time         Lambda
500000     1000.00000        0.00000

Writing checkpoint, step 500000 at Wed Aug 29 01:54:51 2012

Energies (kJ/mol)
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      Potential
5.72530e+03   -2.52945e+02   -3.72877e+04   -4.50918e+03   -3.63246e+04
Kinetic En.   Total Energy  Conserved En.    Temperature Pres. DC (bar)
6.59347e+03   -2.97311e+04   -1.25805e+04    2.93870e+02   -1.54019e+02
Pressure (bar)
-1.82664e+02
<======  ###############  >
<==  A V E R A G E S  ====>
<==  ###############  ======>
Statistics over 500001 steps using 50001 frames
Energies (kJ/mol)
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      Potential
5.74509e+03   -2.52945e+02   -3.72807e+04   -4.52328e+03   -3.63118e+04
Kinetic En.   Total Energy  Conserved En.    Temperature Pres. DC (bar)
6.58722e+03   -2.97246e+04   -1.32539e+04    2.93591e+02   -1.54019e+02
Pressure (bar)
-4.77602e+01

This looks fine.

Mark

; VARIOUS PREPROCESSING OPTIONS
; Preprocessor information: use cpp syntax.
; e.g.: -I/home/joe/doe -I/home/mary/roe
include =
; e.g.: -DI_Want_Cookies -DMe_Too
;define = -DPOSRES

; RUN CONTROL PARAMETERS
integrator = md-vv
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 500000
; For exact run continuation or redoing part of a run
init_step = 0
; Part index is updated automatically on checkpointing (keeps files separate)
simulation_part = 1
; mode for center of mass motion removal
; energy calculation and T/P-coupling frequency
comm-mode = linear
; number of steps for center of mass motion removal
nstcomm = 10
; group(s) for center of mass motion removal
comm-grps =

; LANGEVIN DYNAMICS OPTIONS
; Friction coefficient (amu/ps) and random seed
bd-fric = 91
ld-seed = 1993

; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 100
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10

; TEST PARTICLE INSERTION OPTIONS
rtpi = 0.05

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 0
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file
nstlog = 500
nstenergy = 500
;nstcalcenergy = 1
; Output frequency and precision for xtc file
nstxtcout = 0
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps = System
; Selection of energy groups
energygrps =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 10
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic_molecules = no
; nblist cut-off
rlist = 1.0
; long-range cut-off for switched potentials

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 1.0
; Relative dielectric constant for the medium and the reaction field
epsilon_r = 1
epsilon_rf = 1
; Method for doing Van der Waals
vdwtype = switch
; cut-off lengths
rvdw_switch = 0.8
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = enerpres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Seperate tables between energy group pairs
energygrp_table =
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.10
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-6
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no

; IMPLICIT SOLVENT ALGORITHM
implicit_solvent = no

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
Tcoupl = v-rescale
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau-t = 1
ref-t = xxx
; Pressure coupling
Pcoupl = no
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1
compressibility = 4.5e-5
ref-p = 1.0
; Scaling of reference coordinates, No, All or COM
refcoord_scaling = No
; Random seed for Andersen thermostat
andersen_seed = 815131

; OPTIONS FOR QMMM calculations
QMMM = no
; Groups treated Quantum Mechanically
QMMM-grps =
; QM method
QMmethod =
; QMMM scheme
QMMMscheme = normal
; QM basisset
QMbasis =
; QM charge
QMcharge =
; QM multiplicity
QMmult =
; Surface Hopping
SH =
; CAS space options
CASorbitals =
CASelectrons =
SAon =
SAoff =
SAsteps =
; Scale factor for MM charges
MMChargeScaleFactor = 1
; Optimization of QM subsystem
bOPT =
bTS =

; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing =
; Number of time points to use for specifying annealing in each group
annealing_npoints =
; List of times at the annealing points for each group
annealing_time =
; Temp. at each annealing point, for each group.
annealing_temp =

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = yes
gen-temp = xxx
gen-seed = 173529

; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint_algorithm = lincs
; Do not constrain the start configuration
continuation = yes
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs_iter = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials
morse = no

; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =

; WALLS
; Number of walls, type, atom types, densities and box-z scale factor for Ewald
nwall = 0
wall_type = 9-3
wall_r_linpot = -1
wall_atomtype =
wall_density =
wall_ewald_zfac = 3

; COM PULLING
; Pull type: no, umbrella, constraint or constant_force
pull = no

; NMR refinement stuff
; Distance restraints type: No, Simple or Ensemble
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting = Conservative
; Use sqrt of the time averaged times the instantaneous violation
disre-mixed = no
disre-fc = 1000
disre-tau = 0
; Output frequency for pair distances to energy file
nstdisreout = 100
; Orientation restraints: No or Yes
orire = no
; Orientation restraints force constant and tau for time averaging
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) and S to energy file
nstorireout = 100
; Dihedral angle restraints: No or Yes
dihre = no
dihre-fc = 1000

; Free energy control stuff
free-energy = no
init-lambda = 0
delta-lambda = 0
foreign_lambda =
sc-alpha = 0
sc-power = 0
sc-sigma = 0.3
couple-moltype =
couple-lambda0 = vdw-q
couple-lambda1 = vdw-q
couple-intramol = no

; Non-equilibrium MD stuff
acc-grps =
accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
deform =

; Electric fields
; Format is number of terms (int) and for all terms an amplitude (real)
; and a phase angle (real)
E-x =
E-xt =
E-y =
E-yt =
E-z =
E-zt =

; User defined thingies
user1-grps =
user2-grps =
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0

spc900_equil_NVT.tpr (79.2 KB) spc900_equil_NVT.tpr Mark Abraham, 09/04/2012 06:54 AM
spc900_equil_NVT.tgz (587 KB) spc900_equil_NVT.tgz Input and output files Mark Abraham, 09/04/2012 06:54 AM

Associated revisions

Revision 5f1fd720 (diff)
Added by Mark Abraham almost 5 years ago

Velocity-verlet integrators and v-rescale are broken

Added fatal error and suggestion to upgrade GROMACS version.

Fixes #1000

Change-Id: I22653fafdfd0fb2f1fa772fbd95bc55fcf072d8c

Revision b77fa706 (diff)
Added by Michael Shirts over 4 years ago

Some changes for md-vv extracted from 4.6

a. Fixes for the pressure in MTTK with constraints + dispersion + rerun
  • Dispersion is correctly added in rerun
  • COM motion is removed only on the second half of the timestep.
  • Now can do md-vv + rerun with multiple threads.
  • Now gives exact kinetic energy reruns for everything except MTTK, where the iterative algorithm
    makes exact kinetic energy impossible when nstpcouple == 1.

b. md-vv works with v-rescale and berendsen

c. Fixes a bug when pressure control in md-vv when nstcalcenergy is not a
multiple of nstpcouple or nsttcouple. This bug results in boxes slowly
expanding to unphysical sizes because the virial is
neglected in the second half of the md-vv calculation.

Also discovered that as part of the bug, global energies were being communicated
where they did not need to be when nstpcouple and nsttcouple are > 1 in the case
of md-vv, so redid some of the iteration counting and global communication to fix
this all together. In the process, this simplified some of the iteration counting.

Should fix bugs #1116, #1012, #1000, #1129 in redmine.

Change-Id: I1b628d03ab588c29fef2b8789e61254da49c2b6f

History

#1 Updated by Mark Abraham about 5 years ago

I looked at this again with a system with a single particle in a harmonic potential well, and found the same kind of problem - in release-4-5-patches, md with v-rescale was conservative, and md-vv was not conservative and the sampling was obviously skewed.

#2 Updated by Michael Shirts about 5 years ago

Yes, as I mentioned, this is fixed in 4.6, and I just need to find the time to isolate the fix and pull in from release-4-6. There is very little time during the week for me.

#3 Updated by Mark Abraham almost 5 years ago

  • Status changed from New to In Progress

https://gerrit.gromacs.org/2047 introduces a fatal error to "fix" this

#4 Updated by Michael Shirts almost 5 years ago

Quick thought on this -- is Berendsen already blocked by the 4.5 code? The same problem (and fixes in 4.6) should exist in both cases.

#5 Updated by Mark Abraham almost 5 years ago

Michael Shirts wrote:

Quick thought on this -- is Berendsen already blocked by the 4.5 code? The same problem (and fixes in 4.6) should exist in both cases.

Yes, REMD has blocked Berendsen for some time. My observations in August were that md-vv & v-rescale worked in release-4-6, but clearly that needs further checking now.

#6 Updated by Mark Abraham over 4 years ago

  • Target version changed from 4.5.6 to 4.5.7

Draft https://gerrit.gromacs.org/#/c/2101/ addresses this properly

#7 Updated by Mark Abraham over 4 years ago

  • Status changed from In Progress to Resolved
  • Affected version set to 4.5.5

#8 Updated by Rossen Apostolov almost 4 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF