## Bug #1000

### Possible bug with v-rescale thermostat and md-vv

**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

### Associated revisions

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 over 7 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 over 7 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 about 7 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 about 7 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 about 7 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 about 7 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 almost 7 years ago

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

#### #8 Updated by Rossen Apostolov over 6 years ago

**Status**changed from*Resolved*to*Closed*

Velocity-verlet integrators and v-rescale are broken

Added fatal error and suggestion to upgrade GROMACS version.

Fixes #1000

Change-Id: I22653fafdfd0fb2f1fa772fbd95bc55fcf072d8c