~"user bug report" ## 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 |