Changes of default `nstpcouple` and `nsttcouple` from version 2022 to 2023

User discussions
1

GROMACS version: 2023
GROMACS modification: No

Dear all,

This is not a question, just something I thought may be worth posting for information.

The new 2023 Gromacs made a change in the default nstpcouple and nsttcouple from 10 (v. 2022) to 100 (v. 2023):

Increased default T- and P-coupling intervals
The default maximum values temperature and pressure coupling intervals have been increased from 10 to 100 steps. These values are used when the default value of -1 is specified in the mdp file and a lower value is used when required for accurate integration. The improves the performance of both GPU runs and parallel runs.

Since our mdp was set to defaults, this change had an effect on our simulations, making them unstable.

We learned it accidentally when rerunning a simulation (same structure, topology, mdp) that we have previously run without any issues on v. 2022.2, but we got the following error and, obviously, a wobbly box.

Step 14700  Warning: pressure scaling more than 1%, mu: 0.981371 0.981371 0.99646
Step 14800  Warning: pressure scaling more than 1%, mu: 1.02037 1.02037 1.0013
Step 14900  Warning: pressure scaling more than 1%, mu: 0.979776 0.979776 0.99901

Changing the .mdp input to explicitly declare the T- & P- coupling to be 10, solved the issue.

While this is trivial, I thought it worth mentioning in case anyone else is not reading the changes log when updating to a new version;)

BW,
V

PS. I have to say that there is a noticeable improvement in performance with v. 2023 over 2022 - well done developers:)

2

Thank you for providing this feedback.

I made this change with the idea that this should never cause instabilities when the parameters are set correctly and the system is stable. What temperature and pressure coupling algorithms are you using and what coupling times and compressibility?

There has been a post before about a similar issue, but IIRC there the issue was that the actually compressibility was significantly lower than the value used in the mdp file (which was the standard water value). That caused the period of the Parrinello-Rahman barostat to be much shorter than tau_p, which in turn required a smaller nstpcouple.

3

Hi Berk!

I have not used anything extraordinary. Here is the relevant part of my original MDP file

; Temperature coupling is on
tcoupl          = nose-hoover  			; or V-rescale - modified Berendsen thermostat
tc-grps         = System                ;
tau_t           = 0.1                   ; time constant, in ps
ref_t           = 300                   ; reference temperature, one for each group, in K


; Pressure coupling 
Pcoupl                   = C-rescale	;
Pcoupltype               = semiisotropic 	; separate coupling in xy from z
tau_P                    = 1.0          	; coupling time constant (ps)
compressibility          = 4.5e-5 4.5e-5 	; compressibility (1/bar)
ref_P                    = 1.0 1.0 		; reference P (bar)

We literally were rerunning exactly same system that we have already run on 2022.3 without any issues.

The input parameters from the 2022.3 .log are then:

   tcoupl                         = Nose-Hoover
   nsttcouple                     = 10
   nh-chain-length                = 1
   print-nose-hoover-chain-variables = false
   pcoupl                         = C-rescale
   pcoupltype                     = Semiisotropic
   nstpcouple                     = 10
   tau-p                          = 1
   compressibility (3x3):
      compressibility[    0]={ 4.50000e-05,  0.00000e+00,  0.00000e+00}
      compressibility[    1]={ 0.00000e+00,  4.50000e-05,  0.00000e+00}
      compressibility[    2]={ 0.00000e+00,  0.00000e+00,  4.50000e-05}
   ref-p (3x3):
      ref-p[    0]={ 1.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    1]={ 0.00000e+00,  1.00000e+00,  0.00000e+00}
      ref-p[    2]={ 0.00000e+00,  0.00000e+00,  1.00000e+00}
   refcoord-scaling               = No
   posres-com (3):

and from 2023 are:

   ensemble-temperature-setting   = constant
   ensemble-temperature           = 300
   tcoupl                         = Nose-Hoover
   nsttcouple                     = 5
   nh-chain-length                = 1
   print-nose-hoover-chain-variables = false
   pcoupl                         = C-rescale
   pcoupltype                     = Semiisotropic
   nstpcouple                     = 100
   tau-p                          = 1
   compressibility (3x3):
      compressibility[    0]={ 4.50000e-05,  0.00000e+00,  0.00000e+00}
      compressibility[    1]={ 0.00000e+00,  4.50000e-05,  0.00000e+00}
      compressibility[    2]={ 0.00000e+00,  0.00000e+00,  4.50000e-05}
   ref-p (3x3):
      ref-p[    0]={ 1.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    1]={ 0.00000e+00,  1.00000e+00,  0.00000e+00}
      ref-p[    2]={ 0.00000e+00,  0.00000e+00,  1.00000e+00}

Then I manually add couplings to be 10 into my mdp, the P issue is resolved.

   nsttcouple                     = 10
   nstpcouple                     = 10

Note, my system does need to compress a little on the z-axis. Obviously, with the settings to 10 it was OK, but 100 may be a bit too ambitious. I think you are right that 100 is OK for an equilibrated system (have not tried, as my goal was just to rerun the previous run).

I now know about this change and will make sure to run first a separate equilibration with coupling set to 10, then a separate production run with the larger default coupling.

BW
V

4

A period of 0.1 with Nose-Hoover is too short. You should use at least tau-t=1 ps.

And I would suggest to use tau-p=5 ps. Having the box change on a 1 ps time scale with semi-isotropic could cause issues.

PS: I think your nsttcouple should have been automatically set to 5, not to 100. So this would mean you now increased, not decreased, it to 10.

5

Hi, yes, the tau-t does seem low. My guess is that the reasoning behind it was to set tau-t to be less than tau-p, which is set to the default 1 ps. Yep, increasing both tau-t and tau-p would indeed be better indeed. I have seen and tested the effect of tau-t tau-p on the semiisotropic systems, especially more rigid than the current simulation, and the period had to be increased. Nevertheless, while the set up is not ideal, it does not seem to generate an issue with this current system when run with v. 2022.

I am not sure what is the effect of temperature on the pressure scaling. But changes from nsttcouple=5 (set automatically) to 10 (default of 2022) should not have as dramatic effect as the nstpcouple change from 100 to 10.

6

Just adding for info.

Tried another rerun of a different system that we also did previously with v 2022.3 w/o P scaling issue.
This is just isotropic system.

The .log parameters read:

   ensemble-temperature           = 300
   tcoupl                         = Nose-Hoover
   nsttcouple                     = 25
   nh-chain-length                = 1
   print-nose-hoover-chain-variables = false
   pcoupl                         = Parrinello-Rahman
   pcoupltype                     = Isotropic
   nstpcouple                     = 100
   tau-p                          = 5
   compressibility (3x3):
      compressibility[    0]={ 4.50000e-05,  0.00000e+00,  0.00000e+00}
      compressibility[    1]={ 0.00000e+00,  4.50000e-05,  0.00000e+00}
      compressibility[    2]={ 0.00000e+00,  0.00000e+00,  4.50000e-05}
   ref-p (3x3):
      ref-p[    0]={ 1.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    1]={ 0.00000e+00,  1.00000e+00,  0.00000e+00}
      ref-p[    2]={ 0.00000e+00,  0.00000e+00,  1.00000e+00}
   refcoord-scaling               = No

and we get the Pressure scaling more than 1%.

7

Is this with an equilibrated system?

And what is tau_t in this case?

8

Hi,

I think tau_t was 0.5 ps, not sure why I cannot see it in the .log (?)

Yes, the system could be better equilibrated indeed. I guess, the only point is that the same system with the same mdp run smoothly on v 2022 (and we then only use last parts of the trajectory for any analysis), so we did not blink an eyelid. While in the new version, we definitely need to split into the equilibration part, where we are more careful to declare the constants and then follow with another run with larger/default coupling.

9

That makes sense. I carefully considered that the new defaults should be correct for all production run circumstances. But they could indeed lead to issues under some equilibration conditions.

Now that we have the c-rescale barostat, we advice using the v-rescale thermostat and c-rescale barostat. This avoids all the issues with oscillations of second order coupling algorithms like Nose-Hoover and Parrinello-Rahman.

10

Yep, I guess all I wanted to do is to highlight the change to anyone being lazy like me and going for a long equilibration-into-production run that worked in v 2022, then getting surprised it is not stable.

The barostat: We did this quick check of this same system with C-rescale too - the same issue with pressure for 2023 (I lost the .log file). I am not sure why I did P-R straight away, from the file name it looks to have been an even older .mdp file that would come after equilibration with Berendsen.

11

But this is important for us to know. If a significant fraction of the users has gotten issues with equilibration with these new settings, we should do something about it.

1 Like
12

Hi Valentina and Berk,
We just had a similar issue when simulating a lipid bilayer, using the CHARMM-GUI generated mdp files. For GROMACS, the CHARMM-GUI input generator uses Nose-Hoover thermostat (tau-t=1) and Parrinello-Rahman barostat (tau-p=5), but the nsttcouple and nstpcouple parameters are not explicitly set and, therefore, the default is used (100 steps in 2023). Similarly to Valentina’s simulations, this led to large box fluctuations and instabilities, which solved when we reduced to nstcouple=20 or when we changed to v-rescale and c-rescale. Other users employing the CHARMM-GUI mdp files may experience similar issues.

13

Is this during equilibration or production?

Also semiisotropic pressure coupling will be more sensitive.

14

This is the case for both equilibration and production, as in these two stages CHARMM-GUI mdps do not contain explicit definitions of nsttcouple and nstpcouple.

15

OK, but then I assume the issue shows up with equilibration. We in general do not advise to use Nose-Hoover and Parrinello-Rahman during equilibration.

16

Perhaps we should improve documentation to note that defaults aim to cater for equilibrated systems and equilibration protocols should be done with specific simulation settings suitable for this stem (including type of coupling, tau, etc)?

17

Parrinello-Rahman should not be used during equilibration. grompp warns about this, but can only detect equilibration by checking for generating velocities.

So the issue here is that charmm-gui sets Parrinello-Rahman for equilibration.

I noticed a few outdated details in our documenation about barostats. If more can be improved, please tell me.

1 Like
18

Yes, it would be good to have this in the documentation. But, am I correct that when running equilibration still the automatic coupling will default to 100, independently of the thermo/barostat? Which, even with Berendsen, is too big. So, then you rely on a user to manually set this parameter. I have to honestly admit, while I have been using Gromacs since version 3 (yes I am old), often for materials, this was the first time I added these into the .mdp file.

19

With Berendsen, which you should not use anymore, and c-rescale, the default nstpcouple=min(100, tau_p/(5dt))
With PR it is min(100, tau_p/(25dt))

I would advise to use a tau_p of at least 5 ps, which I suggested to have as default value from the 2024 release onward. I would not think you would get issues with c-rescale and tau_p = 5 ps, but maybe with tau_p = 1 ps you do.

20

Thank you very much for your feedback.

Sorry for one confusing aspect I did not mention: equilibration steps were carried out with Berendsen and production runs with Parrinello-Rahman, both with the default nstpcouple. This is the CHARMM-GUI recommendation. But as 100 steps are used in GMX2023, both stages lead to instabilities.

In another test, I took a 500-ns C-rescale well-equilibrated conformation as starting configuration. I changed to Parrinello-Rahman and left nstpcouple=default (meaning 100 steps). In such case, instabilities are also observed and the simulation crashes after a couple of ns. Thus, instabilities seem to appear for Parrinello-Rahman when using nstpcouple=100, regardless of whether the starting conformation is well equilibrated or not.