Box distortion in simulation of an ion channel with electric field

GROMACS version: 2020.6
GROMACS modification: No

Dear all,

I am currently trying to set up a simulation of an ion channel with an applied electric field of 60 mV for direction z. I now encounter the problem that the simulation box keeps shrinking in the z-axis (see snaphot), resulting in a very flat and distorted system. However, GMX does not output any error messages or warnings during the run. Our group also encountered this problem with Gromacs 2018 on a different ion channel. Runs without electric field are stable and do not shrink.

Can you recommend any way to handle this problem? I read in the forum that a solution might be to switch to isotropic pressure coupling, but since the system contains a membrane, I am not sure about the validity of this approach. Might a longer equilibration help?

Thanks in advance for your advice!
All the best,
Theres

Here a bit of information about the system:
equilibration: 1 ns nvt, 1 ns npt
electric field: 60 mV in z direction
mdrun: 50 ns, 4 runs - in all runs the box shrinks in z-dimension during the simulation

mdp file:

; Run parameters
integrator = md
nsteps = 25000000
dt = 0.002
; Output control
nstxout-compressed = 10000
nstenergy = 10000
nstlog = 10000
; Bond parameters
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid
nstlist = 10
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; Temperature coupling is on
tcoupl = V-rescale
tc-grps = Protein_MOL_POPC Water_Ion_ETR
tau_t = 0.5 0.5
ref_t = 310 310
; Pressure coupling is on
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 2.0
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
; Periodic boundary conditions
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = no
; COM motion removal
comm-mode = Linear
comm-grps = Protein_MOL_POPC Water_Ion_ETR
;Electric fieldselectric field
Electric-field-z = 0.06 0 0 0 ; electric field 0.06 Volt per nm in z direction

This appears to be a rather obvious bug. Please try isotropic pressure coupling to see if you at least get sensible coordinates. If you do, this points to a problem specifically in semiisotropic coupling, which is important information for correcting the problem.

Thanks a lot for the swift answer!

I have now rerun the simulation with isotropic pressure coupling. Unfortunately, the simulation does not look a lot better at the end.

As already mentioned, the system is perfectly stable in the free MD simulation. Can you recommend any next steps (perhaps check the lipid density of the membrane or equilibrate longer?) Or should I report my problem

By the way, I forgot to mention that I am using amber99sb in combination with Berger lipids.

Thank you again very much for your help!

Please file a bug report on the GROMACS GitLab site. The fact that the systems distort and compress so badly without failing is suspicious.

I will do that, thank you for the support!