GROMACS version: 5.0.2 and 2021.2
GROMACS modification: No
I am studying adsorbtion of organic surfactants onto the surface of inorganic crystal using replica-exchange simulations in NVT ensemble (replicas between 300 and 2200K). All interactions are treated classically, and the crystal is modelled using only nonbonded interactions (Coulomb and LJ). Everything was working well in GROMACS 5.0.2, but now I am switching to a newer 2021.2 version and there are some issues. Simulations in the highest temperature replicas blow up giving LINCS warnings. I tried to run a separate plain MD simulation at the highest temperature and got the same issue.
I’ve compared mdout.mdp files generated by the two versions of Gromacs. There are no differences that can cause the problem. I am using altered masses, so my frist thought was that the new version of grompp takes masses from the [ atomtypes ] section instead of [ atoms ]. But it appears to be not the cause of the problem.
Another strange observation is that when visualizing the new trajectory at 2200K in VMD I see atoms that slightly move out of the simulation box for some time. This is very suspicious to me, since no periodic gathering was applied. In case of older version everything was fine, all the atoms stayed within the simulation box. Can someone comment on this?
What changes have been made since version 5.0.2 that may alter the stability of MD simulations and can cause the observed problem?
Thank you,
Andriy
I have checked tpr files that were generated by the two versions of grompp. They are identical, except for random seed and buffer size (- denotes GROMACS 5.0.2, and + denotes GROMACS 2021.4):
$ diff -u md_tpr.dat …/…/REMD1-test2-3/t119/md_tpr.dat
— md_tpr.dat 2022-03-18 14:24:15.000000000 +0100
+++ …/…/REMD1-test2-3/t119/md_tpr.dat 2022-03-18 14:22:01.000000000 +0100
@@ -10,7 +10,7 @@
comm-mode = None
nstcomm = 0
bd-fric = 0
- ld-seed = 1591947263
emtol = 10
emstep = 0.01
niter = 20
@@ -176,7 +176,7 @@
bF = not present
natoms = 10757
lambda = 0.000000e+00
- buffer size = 388782
topology:
name=“FAPbBr3_NMe3C2PC5_solvated”
#atoms = 10757
Therefore, it seems that the problem is not with the parameters or topology files, but there is some bug in mdrun, which shows up in my simulations.
In 2019, the concept of update groups was introduced, which may potentially affect the stability of my MD simulations. I tried to disable it by setting the corresponding environment variable GMX_NO_UPDATEGROUPS to 1. The problem, however, persists.
Here is the content of my mdp-file:
title = FAPbBr3_NMe3C2PC5_in_toluene
define = -DPOSRES_REMD
constraints = h-bonds
constraint-algorithm = LINCS
integrator = md
dt = 0.001 ; ps !
nsteps = 10000000 ; 10 ns
nstxout = 20000 ; Save coordinates
nstlog = 20000
nstenergy = 20000
nstvout = 20000
nstfout = 0
; mode for center of mass motion removal
comm-mode = None
;
xtc-precision = 1000
xtc-grps = Brom Lead FA LIG1
nstxtcout = 4000
;
cutoff-scheme = Verlet
nstlist = 10 ; frequency to update neighbors list
ns_type = grid ; method to search for neighbors (cell)
rlist = 1.00 ; cut-off distance for the neighbors’ list
; Method for doing Van der Waals =
vdw-type = cut-off
; cut-off lengths =
rvdw-switch = 0.0
rvdw = 1.00 ; cut-off for vdW
; Method for doing electrostatics =
coulombtype = PME; user; cut-off; Reaction-Field; Generalized-Reaction-Field ;
rcoulomb-switch = 0
rcoulomb = 1.00 ; cut-off for Coulomb force
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon-r = 1
; Spacing for the PME/PPPM FFT grid =
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used =
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters =
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; Temperature coupling
Tcoupl = nose-hoover
tc-grps = system
tau_t = 0.5 ; relx. time of T. Real time is ~50 times longer. 0.05 is probably OK
; ref_t = 2000.000 ; Temperature defined in the end of the file
; Pressure coupling is
pcoupl = no ; Parrinello-Rahman ; Berendsen
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 0 1.1e-4 ;box is fixed in x,y directions
ref_p = 1.0 1.0
refcoord-scaling = no ; com ;do not scale origins of each restrain but the center of their mass
; Generate velocites?
gen_vel = no ;yes
gen_temp = 2000.000
gen_seed = 173529
ref_t = 2200.00
Besides the usual force-field terms, my simulation uses harmonic restraints, and flat-bottomed restraints with spherical and layer geometries.
Andriy