GROMACS version: 2025.3
GROMACS modification: No
I am trying to create a box filled with ionic liquid as solvent and some simple lipids. Force field is OPLS-AA. I tried this with lipids and water (TIP4P) first, and it worked fine. However, although my ionic liquid/lipid mix is stable through energy minimisation and NVT/NPT, there are huge void spaces in the liquid which will not go away during NPT. This prevents me from being able to run MD and see the motion of lipids through the IL, as the box just explodes from negative pressure.
I have tried Berensden, C-rescale, and Parrinello-Rahman barostats. Isotropic/semiisotropic made no difference. I have also tried increasing the reference pressure to “squish” the box, decreasing tau_p, increasing dt (in case the system just needs more time to relax), “annealing” the system by increasing then decreasing temperature in NVT, and combinations of the above.
The box is slowly shrinking with decreasing tau_p, but it is taking days of runtime and the voids arent actually disappearing.
The box was initially packed with Packmol, and I had to increase the box size from my calculated density, because the system was failing energy minimisation. EM, NVT, and current NPT setups are below.
____________________________________________________
EM (3 stage, this was the final):
integrator = steep ; can keep steep or switch to conjugate gradient
emtol = 500.0 ; lower than before
emstep = 0.001 ; smaller steps
nsteps = 50000 ; or more if needed
nstlist = 1
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.2
rvdw = 1.2
pbc = xyz
NVT (non-annealing):
; NVT equilibration - match paper exactly
; Velocity-rescaling weak-coupling, 300 K, tau_T = 0.1 ps
; Paper: Temperature control via velocity-rescaling weak-coupling, 300 K
title = NVT equilibration
; Run parameters
integrator = md
dt = 0.002
nsteps = 50000
nstxout = 500
nstvout = 500
nstenergy = 500
nstlog = 500
; Output control
energygrps = System
; Neighbor list
cutoff-scheme = Verlet
nstlist = 10
pbc = xyz
; Electrostatics - PME
coulombtype = PME
rcoulomb = 1.0
; Van der Waals - cutoff per force field
vdwtype = Cut-off
rvdw = 1.0
; Temperature coupling - velocity-rescaling weak-coupling
; Paper: tau_T = 0.1 ps, 300 K
tcoupl = V-rescale
tc-grps = System
tau_t = 0.1
ref_t = 300
; Pressure - no barostat for NVT
pcoupl = no
; Constraints - LINCS for bonds, SETTLE for water
; Paper: LINCS for bond lengths, SHAKE for H-bonds, SETTLE for water
; GROMACS uses LINCS for h-bonds; use constraint-algorithm = shake to match paper
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
NPT (most recent version, tightest tau_p):
; NPT - Based on glycine_hso4/paper_replication protocols
; tau_p 1 ps
title = NPT with head group restraints
define = -DPOSRES_HEAD
; Run parameters - 200 ps (glycine paper uses 2 ns; extend if needed)
integrator = md
dt = 0.002
nsteps = 1000000
nstxout = 500
nstvout = 500
nstenergy = 500
nstlog = 500
; Output control
energygrps = System
; Neighbor list (paper: nstlist 20)
cutoff-scheme = Verlet
nstlist = 20
pbc = xyz
; Electrostatics - PME (paper: rcoulomb 1.2, ewald-rtol 1e-6)
coulombtype = PME
rcoulomb = 1.2
fourierspacing = 0.12
ewald-rtol = 1e-6
; Van der Waals (paper: rvdw 1.2)
vdwtype = Cut-off
rvdw = 1.2
; Temperature coupling (paper: tau_t 1.0)
tcoupl = V-rescale
tc-grps = System
tau_t = 1.0
ref_t = 300
; Pressure coupling - glycine_hso4 C-rescale protocol
pcoupl = C-rescale
pcoupltype = isotropic
tau_p = 0.75
ref_p = 1.0
compressibility = 4.5e-5
refcoord_scaling = com
; Constraints (v4 settings for stability)
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 6
lincs-iter = 2
lincs-warnangle = 30
______________________________________________________
The simulation box looks like this
The “branches” of liquid are very stable, and stay fairly rigid over time.
Is this a known problem for viscous, charged liquids, and is there a more straightforward way to overcome this void space?
Thanks