Z‑dimension of the box decreases

GROMACS version:2018
GROMACS modification: Yes/No
Here post your question
“Dear all, I am using the Martini 3‑IDP force field to simulate a system consisting of a spherical protein condensate and a cell membrane (DIPC/DIPS) with initial box dimensions of 30 x 30 x 30 nm, built using the command insane -f lk.gro -pbc rectangular -salt 0.15 -x 30 -y 30 -z 30 -sol W -o system.gro -u DIPC:50 -u DIPS:50 -l DIPC:50 -l DIPS:50 -center -dm 9.0. During the production MD run, I observe that the Z‑dimension of the box continuously decreases while the XY area increases, despite using Pcoupltype = semiisotropic. I have already tried changing tau_p from 12.0 to 5.0 and compressibility from 3e‑4 3e‑4 to 4.5e‑5 4.5e‑5, but the problem persists. Could anyone suggest how to resolve this box compression issue? Thank you.”
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.

title = Martini

; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.

integrator = md
dt = 0.01
nsteps = 1000000000
nstcomm = 100
comm-grps =

nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxout-compressed = 10000
compressed-x-precision = 300
compressed-x-grps =
energygrps = system

; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.

cutoff-scheme = Verlet
nstlist = 30
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005

coulombtype = reaction-field
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable water)
epsilon_rf = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1

; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.

tcoupl = V-rescale
tc-grps = Protein Bilayer SOL
tau_t = 1.0 1.0 1.0
ref_t = 295 295 295
Pcoupl = Parrinello-rahman
Pcoupltype = semiisotropic
tau_p = 12 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0

gen_vel = no
gen_temp = 295
gen_seed = -1

; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.

constraints = none
constraint_algorithm = Lincs

That is not “despite” semi-isotropic coupling, as that is exactly what allows this. One reason could be that your protein is not well embedded in the membrane, causing large lateral forces.

tau_p=1 is far too short for Parrinello-Rahman. And PR-coupling can be unstable, especially during equilibration. We recommend using C-rescale with tau_p at least 5 ps.

Thank you for your reply. Unfortunately, we uniformly use GROMACS 2018, which does not support C‑rescale, and my tau‑p is set to 12. However, I have now resolved this issue. The reason was that the membrane built by insane is initially too regular and dense, causing the membrane area to expand during the simulation. As a result, the box area in the xy plane also increases, while the box length along the z axis decreases accordingly. The solution was to first simulate the membrane alone in a solvent environment until the membrane had expanded to its maximum extent, and then place the protein and the membrane into a new box, which solved the problem. Finally, the initial state of my system is as follows: