Results of MTTK barostat diffs from other barostats
Summary
We found the density obtained with the MTTK barostat differs from that obtained with other barostats in GROMACS in 2024.1.
GROMACS version
2024.1
Steps to reproduce
The attachment in GoogleDrive contains input files used for our testing. Except for the files labeled as 2018.8, all others are input files for version 2024.1. They include topology files, initial configurations, parameter files, and scripts submitted to the cluster.
https://drive.google.com/file/d/1ZpfznDkV1dtxiIbGMlje_yKofWevwDvG/view?usp=sharing
Current bug behavior
In GROMACS version 2018.8, the density obtained using all barostats, including MTTK, is approximately 1.018 g/cm³. However, in the latest version 2024.1, the density obtained with the MTTK barostat is 1.043 g/cm³.
Density | |
---|---|
GROMACS2018.8 Velocity-Verlet MTTK | 1.01785±4.29353E-4 |
GROMACS2024.1 Velocity-Verlet MTTK legacy | 1.04337±4.32314E-4 |
GROMACS2024.1 Velocity-Verlet MTTK modular | 1.04318±8.92863E-5 |
GROMACS2024.1 Velocity-Verlet C-rescale legacy | 1.01831±3.45471E-4 |
GROMACS2024.1 Velocity-Verlet C-rescale modular | 1.01788±3.5554E-4 |
GROMACS2024.1 Leap-frog C-rescale legacy | 1.01731±2.74053E-4 |
GROMACS2024.1 Leap-frog C-rescale modular | 1.01789±5.45831E-4 |
GROMACS2024.1 Leap-frog PR legacy | 1.01788±2.6744E-4 |
GROMACS2024.1 Leap-frog PR modular | 1.01771±5.76678E-4 |
We tested various barostats, including C-rescale, MTTK, and Parrinello-Rahman, under both leap-frog and velocity-verlet algorithms in both the legacy and modular simulators. It can be observed inconsistent density results are obtained only when using the MTTK pressure controller.
Subsequently, we tested multiple historical versions of GROMACS between 2018.8 and 2024.1. Until version 2020.7, the results with the MTTK pressure controller remained consistent with other barostats. Unfortunately, we were unable to successfully compile any versions from 2021. However, starting from version 2022, inconsistent results were obtained with the MTTK barostat.
For developers: Why is this important?
The converged values of density should be the same for different barostats that was proved to exactly reproduce the NPT ensembles. The inconsistency of different methods in GROMACS is confusing for users.
We look forward to the development team carefully analyzing our feedback and promptly addressing this issue. We appreciate their support and attention to our concerns.
input mdp example
title = MET NPT
define = -DFLEXIBLE
integrator = md-vv
nsteps = 50000000
dt = 0.0002
continuation = no
nstxout = 1000
nstvout = 0
nstenergy = 200
nstlog = 200
constraint_algorithm = lincs
constraints = none
cutoff-scheme = Verlet
ns_type = grid
nstlist = 10
rlist = 0.7
rcoulomb = 0.9
rvdw = 0.7
coulombtype = PME
fourierspacing = 0.16
tcoupl = nose-hoover
tc-grps = System
nh-chain-length = 4
tau_t = 2
ref_t = 298.15
pcoupl = MTTK
pcoupltype = isotropic
tau_p = 4.0
ref_p = 1.013
compressibility = 4.5e-5
refcoord-scaling = all
pbc = xyz
DispCorr = EnerPres
gen_vel = yes
nsttcouple = 1
nstpcouple = 1
Update
Furthermore, we have tested the spc/fw model. The simulation of the spc/fw model in the literature yielded a density of 1.012 (±0.016) g/cm³.
Here is the literature link: https://doi.org/10.1063/1.2136877.
The result of C-rescale is still diff from MTTK barostat, but more close to the reference value.
Density | |
---|---|
Origin Work | 1.012±0.016 |
GROMACS2024.1 Velocity-Verlet MTTK | 1.03496±1.82012E-4 |
GROMACS2024.1 Velocity-Verlet C-rescale | 1.00938±4.43439E-4 |
GROMACS2024.1 Stochastic Dynamics C-rescale | 1.00911±2.69820E-4 |