Parrinello-Rahman instability in Gromacs 2023

User discussions
1

GROMACS version:2023
GROMACS modification: No

Hello everyone!
Recently i’ve moved to Gromacs 2023 and really happy about it’s performance. But occasionally, it happens to segfault on simulations in N.-H. & P.-R. NPT of one of my systems. It warns with “Pressure scaling more than 1%”, despite system is seems to be well equilibrated in Berendsen NPT, and after 5-6 ns crashes with segfault. I’ve tried to grompp and mdrun from exact same files on other machine with Gromacs 2021.2, and faced no errors and warnings on 10x simulations length. Also, the cell dimensions simulated with older Gromacs seem to be less fluctuant.
Maybe someone had the same problem and knows how to solve it? For now, the only option i have is to downgrade to older versions of Gromacs.

I dont use any of experimental env. variables, run gmx mdrun as it is without specifying anything more than a simulation name (since automatic mdrun performance choice works best on my hardware).

Hardware and System

GPU NVIDIA GeForce 2080Ti
CPU AMD Ryzen 7 3700X
16 Gb RAM
Ubuntu 22.04

Here is my gmx --version output

GROMACS version:    2023
Precision:          mixed
Memory model:       64 bit
MPI library:        thread_mpi
OpenMP support:     enabled (GMX_OPENMP_MAX_THREADS = 128)
GPU support:        CUDA
NB cluster size:    8
SIMD instructions:  AVX2_256
CPU FFT library:    fftw-3.3.8-sse2-avx-avx2-avx2_128
GPU FFT library:    cuFFT
Multi-GPU FFT:      none
RDTSCP usage:       enabled
TNG support:        enabled
Hwloc support:      disabled
Tracing support:    disabled
C compiler:         /usr/bin/cc GNU 11.3.0
C compiler flags:   -fexcess-precision=fast -funroll-all-loops -mavx2 -mfma -Wno-missing-field-initializers -O3 -DNDEBUG
C++ compiler:       /usr/bin/c++ GNU 11.3.0
C++ compiler flags: -fexcess-precision=fast -funroll-all-loops -mavx2 -mfma -Wno-missing-field-initializers -Wno-cast-function-type-strict -fopenmp -O3 -DNDEBUG
BLAS library:       External - detected on the system
LAPACK library:     External - detected on the system
CUDA compiler:      /usr/local/cuda-12.0/bin/nvcc nvcc: NVIDIA (R) Cuda compiler driver;Copyright (c) 2005-2023 NVIDIA Corporation;Built on Fri_Jan__6_16:45:21_PST_2023;Cuda compilation tools, release 12.0, V12.0.140;Build cuda_12.0.r12.0/compiler.32267302_0
CUDA compiler flags:-std=c++17;--generate-code=arch=compute_50,code=sm_50;--generate-code=arch=compute_52,code=sm_52;--generate-code=arch=compute_60,code=sm_60;--generate-code=arch=compute_61,code=sm_61;--generate-code=arch=compute_70,code=sm_70;--generate-code=arch=compute_75,code=sm_75;--generate-code=arch=compute_80,code=sm_80;--generate-code=arch=compute_86,code=sm_86;--generate-code=arch=compute_89,code=sm_89;--generate-code=arch=compute_90,code=sm_90;-Wno-deprecated-gpu-targets;--generate-code=arch=compute_53,code=sm_53;--generate-code=arch=compute_80,code=sm_80;-use_fast_math;-Xptxas;-warn-double-usage;-Xptxas;-Werror;-D_FORCE_INLINES;-fexcess-precision=fast -funroll-all-loops -mavx2 -mfma -Wno-missing-field-initializers -Wno-cast-function-type-strict -fopenmp -O3 -DNDEBUG
CUDA driver:        12.0
CUDA runtime:       12.0

Nvidia Driver versions

NVIDIA-SMI 525.85.12    Driver Version: 525.85.12    CUDA Version: 12.0

mdp file

integrator              = md
dt                      = 0.002
nsteps                  = 500000000         ;1000ns
nstlog                  = 1000
nstxout-compressed      = 10000
nstvout                 = 10000
nstfout                 = 10000
nstcalcenergy           = 100
nstenergy               = 1000
;
cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2
coulombtype             = pme
rcoulomb                = 1.2
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw_switch             = 1.0
rvdw                    = 1.2
;
tcoupl                  = Nose-Hoover
tc_grps                 = MEMB   SOLV
tau_t                   = 1.0	1.0
ref_t                   = 310	310
;
pcoupl                  = Parrinello-Rahman
pcoupltype              = semiisotropic
tau_p                   = 5.0
compressibility         = 4.5e-5  4.5e-5
ref_p                   = 1.0     1.0
;
constraints             = h-bonds
constraint_algorithm    = LINCS
continuation            = yes
;
nstcomm                 = 100
comm_mode               = linear
comm_grps               = MEMB   SOLV
;
refcoord_scaling        = com

Thank you!

2

The only relevant difference I can think of is that we increased the default nsttcouple and nstpcouple values. I think they should end up both at 100 in your run with 2023. You can check this in your log file. They are likely smaller in your 2022 runs (check the log file). Could you put the 2022 values for nsttcouple and nstpcouple in your mdp file for 2023, run grompp and mdrun and check if that solves the issue?

Please report back here so we can follow up on this.

1 Like
3

Note that I would not actually expect issues with the new default values.These will give at least 20 integration steps per period. What could be the issue is that you actual compressibility is much lower than the values (for water) you use for the barostat. This could lead to a (too) short period. What period of volume fluctuations do you observe in the runs with version 2022 (you might need to do a short run with more frequent energy output to estimate the period)?

1 Like
4

Another possibly relevant change in GROMACS 2023 is using the GPU-resident mode by default. You can try adding -update cpu flag to gmx mdrun to revert to the old (force-offload) mode and see if it improves the situation.

1 Like
5

Thank you for your responce!
Lowering nstt(p)couple to 10 solved the issue for Gromacs 2023. Values of nstt(p)couple in crashed runs were 100.
I ran 100 ps simulation by Gromacs 2021.2 as you’ve suggested to and observed 3-5 ps period for volume fluctuations.

gmx energy summary

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Temperature                 310.003       0.02     3.2832  -0.021334  (K)
Volume                      223.317       0.19   0.899733   0.836876  (nm^3)

Temperature dependent fluctuation properties at T = 310.003.

Heat capacities obtained from fluctuations do *not* include
quantum corrections. If you want to get a more accurate estimate
please use the g_dos program.

WARNING: Please verify that your simulations are converged and perform
a block-averaging error analysis (not implemented in g_energy yet)
Volume                                   = 3.20043e-05 m^3/mol
Isothermal Compressibility Kappa         = 8.45753e-10 (m^3/J)
Adiabatic bulk modulus                   = 1.18238e+09 (J/m^3)
6

Good that you got it too work.

But even with a period of 3 ps there are still 15 integration steps per period of the barostat, which I think should be stable.

Is the simulation also stable when you reduce the compressibility values in your mdp file by a factor 1.5.

7

Yes, i’ve tried 3.0e-5 (1.5 times lower) compressibility for both membrane and water and got only one “Pressure scaling” warning in the very beginning (60000 step) of the simulation. If i keep one of the groups (memb or solv) compressibilities at 4.5e-5, system is less, but still instable. Lowering compressibilities to 1.5e-5 leads to no warnings at all, but i doubt about how physical these values are. Of course, simulations were carried out with nst(p/t)couple equal to 100.
Also, i’ve tried to calculate compressibility as given in How to determine compressibility? thread, but came up with numbers twice bigger than 4.5e-5 bar-1. System exploding almost instantaneously with that values being set.
Thank you for your suggestions and overall help.

8

Little Update
On 3e-5 bar-1 system is still exploding, but it takes much more time, since calculation is slowing down because of pair search takes much more time.
On 1.5e-5 bar-1 simulation is stable.

9

Decreasing the compressibility is mathematcally identical to increasing tau_p by the same factor, as these parameter only enter through their ratio in the barostat. The only difference in GROMACS is that nstpcouple is determined based on tau_p.

So the question is still if the actual x/y compressibility of your system is low or the procedure for the determining nstpcouple in GROMACS is not conservative enough.

For determine compressibility in a membrane system you need to look at x/y-plane area fluctuations, not volume fluctuations. The compressiblity along z will be close enough to that of pure water that you can keep that value.

10

Hello!
X/Y plane area fluctuations period is ~5 ps. Membrane normal oriented along Z. You can see plots here.