GROMACS version: 2024.3
GROMACS modification: No
I’m running a GROMACS simulation of a mixed methanol/water solvent system. I added 2000 methanol molecules using insert-molecules and 2500 water molecules using solvate. The initial box density seems reasonable.
However, after the NPT equilibration step, the box size increases and the density drops significantly to around 620 kg/m³, while the expected density of a methanol/water mixture should be above 791 kg/m³.
🔹 What could be causing this drop in density, and how can I fix it? Should I revisit how I insert the solvent molecules, or adjust my NPT equilibration settings (e.g., force field choice, barostat parameters, compressibility)?
My protocol is:
gmx editconf -f TFA.gro -o box.gro -c -box 6.0 6.0 6.0
gmx insert-molecules -f box.gro -ci methanol.gro -nmol 2000 -o box_methanol.gro
gmx solvate -cp box_methanol.gro -cs spc216.gro -o solvated.gro -p topol.top
gmx grompp -f EM.mdp -c solvated.gro -p topol.top -o em.tpr -n index.ndx
gmx mdrun -deffnm em -v
gmx grompp -f NVT.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -n index.ndx
gmx mdrun -deffnm nvt -v
gmx grompp -f NPT.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -n index.ndx
gmx mdrun -deffnm npt -v
NPT.mdp file:
title = TFA in methanol NPT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 500000 = 1000 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = yes ; continuing from NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 2 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = BER Solvent_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = C-rescale ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; velocity generation off after NVT
After NPT step, box size increase from 6x6x6 to 6.8x6.8x6.8 and density decrease to 620 kg/m³
What wrong with my protocol?