What is the recommended value for tau-p and tau-t with GROMACS2023.0 when using PR and NH

User discussions
1

GROMACS version: 2023.0
GROMACS modification: No

Hey there,

After consulting the GROMACS 2023 documentation
and this post on GitLab Increase the default values of nsttcouple and nstpcouple (!2624) · Merge requests · GROMACS / GROMACS · GitLab, I noticed that if the value -1 is used for nsttcouple and nstpcouple, their default value will be adjusted to 100 (if the tau-p and tau-t values allow for this). In the GROMACS versions until 2022 both values were set to 10 if a value of -1 was selected (see GROMACS 2022 documentation.

Additionally, for first-order coupling (such as Berendsen or V-rescale), the tau-p and tau-t values must be at least 5 times greater than (nstpcouple * dt) or (nsttcoupl * dt). For second-order coupling (such as Parrinello-Rahman or Nosé-Hoover), the tau-p and tau-t values must be at least 20 times greater than (nstpcouple * dt) or (nsttcouple * dt). See Strange Fluctuations with NH+PR coupling in GROMACS 2018.3 - Redmine #2749 (#2749) · Issues · GROMACS / GROMACS · GitLab and the GROMACS 2023 documentation.

Furthermore, when using the Parrinello-Rahman barostat and Nosé-Hoover thermostat, the tau-p value must be at least twice the value of tau-t to prevent resonance, as documented here Strange Fluctuations with NH+PR coupling in GROMACS 2018.3 - Redmine #2749 (#2749) · Issues · GROMACS / GROMACS · GitLab.

Based on these guidelines, the minimum tau-t value would be 1 ns for V-rescale and 4 ns for Nosé-Hoover, respectively, as suggested by @hess here Increase the default values of nsttcouple and nstpcouple (!2624) · Merge requests · GROMACS / GROMACS · GitLab. Using the Parrinello-Rahman barostat would result in a minimum tau-p of 8 ns. However, these values seem higher than what I’m used to, and I’m not sure if they could cause my simulation to blow up.

I am simulating a membrane protein (G protein-coupled receptor) with a ligand in a POPC bilayer membrane.

My previous parameters in the .mdp file were:

tcoupl = Nose-Hoover
tc_grps = PROT_LIG MEMB SOLV
tau-t = 1.0 
ref-t = 310 310 310
nsttcouple = 10

pcoupl = Parrinello-Rahman
pcoupletype = semiisotropic
tau-p = 5.0
compressibility = 4.5e-5 4.5e-5
ref-p = 1.0 1.0 
nstcouple = 10

I would change these parameters to:

tcoupl = Nose-Hoover
tc_grps = PROT_LIG MEMB SOLV
tau-t = 4.0 
ref-t = 310 310 310
nsttcouple = 100

pcoupl = Parrinello-Rahman
pcoupletype = semiisotropic
tau-p = 8.0
compressibility = 4.5e-5 4.5e-5
ref-p = 1.0 1.0 
nstcouple = 100

However, I am concerned that the large values for tau-t and tau-p, combined with the Nosé-Hoover and Parrinello-Rahman methods, respectively, could cause the simulation to blow up. Additionally, I am unsure if the tau-p/tau-t ratio is too small for the NH thermostat and PR barostat to avoid resonance. Intuitively, I would increase the tau-p value, but I am worried this would further increase the risk of the simulation failing.

@hess during the tests to increase the default value of nsttcouple and nstpcouple, were the tau-t and tau-p values also tested? Is there any recommendation especially for the combination of NH thermostat and PR barostat?

Thanks
Eddie

2

Firstly, I assume you mean picoseconds and not nanoseconds, right?

Apart from the warning when tau_p < 2 * tau_t, all other comments and changes only pertain to the coupling intervals or time-step, if you prefer to call it that way. Those changes in no way affect the coupling tau values. If your tau values are small, grompp will automatically reduce the nst parameter such the integration accuracy requirements are fulfilled. If you manually set the nst values and they do not obey the integration accuracy criteria, grompp will generate a warning.

So you can (still) set tau_t and tau_p as you like. The risk of the system blowing up is always with too small values of tau when using second order coupling, never with too large values.

I think tau_t=1 ps is on the small side, as this adds temperature oscillations with a period of 1 ps. I would recommend a value of at least 2 ps. But 1 ps often also works fine.

Nowadays the recommendation is to use the v-rescale thermostat and the c-rescale barostat. These have non of the disadvantages of second order coupling and have no disadvantages themselves.

1 Like
3

Thank you very much for your insightful answer!

Yes, I meant picoseconds not nanoseconds.

Perfect than I will try v-rescale as a thermostat and c-rescale as a barostat.

For the combination of v-rescale as thermostat and c-rescale as barostat is there any condition for the ratio of tau_p/tau_t that needs to be fulfilled (such as tau_p >= 2* tau_t when using NH and PR)?

I would really like to use the increased values for nsttcouple and nstpcouple so I would increase my tau_p and tau_p value respectively.

Regarding c-rescale, I could find the option in the in the GROMACS 2023 documentation but I couldn’t find any information about the algorithm in the current reference manual.

I assume though it is a first order coupling.

Based on that, first order couplings should be at least 5 times higher than (nstpcouple *dt) or (nsttcouple * dt) I would assume that with a time step of 2 femtoseconds the tau_p and tau_t should be at least 1 picosecond.

Would it be advisable to use the same value for tau_p as for tau_t? From my experience the tau_p value was always a multiple larger than the tau_t. But to be fair I never knew the explicit reason for that and went with the values generated and recommended by CHARMM-GUI.

Also, I recall reading, that higher values for tau_p and tau_t are preferred because they reduce the “bias” that is introduced by the barostat and thermostat. But again, I do not really know which value for tau_p and tau_t would be considered too high or too low. I basically just follow the recommendations in the GROMACS manual to determine the smallest accepted value.

Is there any guideline or a rule of thumb to determine the “right” value for tau_p and tau_t depending of the corresponding barostat and thermostat, besides the recommendation for the smallest accepted values in the GROMACS manual? I searched all the literature I could find but did not get a satisfying answers.

Based on you input I would adjust my mdp file as follows:


dt = 0.002

tcoupl = v-rescale
tc_grps = PROT_LIG MEMB SOLV
tau-t = 1.0
ref-t = 310 310 310
nsttcouple = -1 # It should be adjusted to 100

pcoupl = c-rescale
pcoupletype = semiisotropic
tau-p = 1.0
compressibility = 4.5e-5 4.5e-5
ref-p = 1.0 1.0
nstcouple = -1 # It should be adjusted to 100

Key question: which tau_p and tau_t values are recommended to run a simulation of a membrane bound protein using v-rescale as thermostat and c-rescale as a barostat using GROMACS 2023.0?

4

C-rescale is described, including reference, in the manual here:
https://manual.gromacs.org/current/reference-manual/algorithms/molecular-dynamics.html#stochastic-cell-rescaling

First order coupling algorithms usually do not have issues with stability due to oscillations, so choosing the value for tau is not so critical. But tau_t needs to be small enough to effectively remove energy drift due to integration errors. tau_p should not be too small, as scaling the volume too fast might cause instabilities.

My advice would be tau_t=1 ps and tau_p=5 ps for nearly all atomistic simulations.

3 Likes
5

Thanks again for your helpful answers and your clear recommendation. Much appreciated!

6

Hi! Very helpful explanation. But I wonder if you have any suggestions on tau_t and tau_p for coarse grained simulations using v-rescale as a thermostat and c-rescale as a barostat. And can I still use compressibility = 4.5e-5 for my 230k-atom system that is polymers solvated with martini3 water?
Many thanks,
Xufan

The whole .mdp file is
title = production run
;;;define = -DPOSRES ; position restrain
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 100000000 ; 20 * 100000000 = 2000 ns
dt = 0.02 ; 20 fs, 30 fs explodes
; Output control
nstxtcout = 50000 ; save coordinates every 500.0 ps
nstvout = 0 ; save velocities every 1.0 ps
nstenergy = 0 ; save energies every 100.0 ps
nstlog = 50000 ; update log file every 500.0 ps
; Bond parameters
continuation = yes ; Restarting after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; none bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; 20 fs, largely irrelevant with Verlet scheme
verlet-buffer-tolerance = 0.005
; nonbonded parameters
rcoulomb = 1.1 ; short-range electrostatic cutoff (in nm)
epsilon_r = 15 ; 2.5 (with polarizable water)
rvdw = 1.1 ; short-range van der Waals cutoff (in nm)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Solute W_Ion ; two coupling groups - more accurate
tau_t = 1.0 1.0 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 1.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no
; com
refcoord_scaling = com
comm-mode = Linear ; Remove center of mass translation
; comm-grps = Protein Non-Protein ; group(s) for center of mass motion removal, default is the whole system

7

I don’t know what the compressibility of Martini models is, ask on the Martini forum/list.

I would think one would want to increase tau_t and tau_p by roughly the factor of the difference in time step. But then one needs to check that the actual temperature is sufficiently close to the reference temperature.