PBC = xy with pressure coupling gives erroneous results

User discussions
, ,
1

GROMACS version: 2024.1
GROMACS modification: No

Good time of the day.

I am trying to simulate a membrane bilayer’s interaction with anti-microbial peptides. Normal PBC condition simulations work just fine, but due to the nature of AMP’s, I’d like to keep them to one side of the membrane. I am trying to use pbc = xy in the mdp parameters, but after a while it gives erroneous coordinates where an atom, or atoms travel way past the simulation box.

In my most recent attempt, this was the error:
“An atom is beyond the wall: coordinates inf inf -inf, distance -1.100000”

This seems to happen only when pressure coupling is enabled.

I’d like to know if I’m doing something wrong, or if there are better ways of achieving my goal?

Thank you very much in advance for any tips you may have.

I am using CG martini lipids and forcefield .Here are the contents of the .mdp file:
title = Martini

integrator = md
dt = 0.02
nsteps = 30000000
nstcomm = 100
comm-grps =

nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 600000 ; Output frequency for energies to log file
nstenergy = 60000 ; Output frequency for energies to energy file
nstxtcout = 600000 ; Output frequency for .xtc file
xtc_precision = 100
xtc-grps = POPC POPG W
energygrps =

cutoff-scheme = Verlet
nstlist = 150
ns_type = grid
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 ;(for use with Verlet-pairlist)
vdw-modifier = Potential-shift-verlet
rvdw = 1.1 ;(for use with Verlet-pairlist)

tcoupl = v-rescale
tc-grps = POPC POPG W
tau_t = 1.0 1.0 1.0
ref_t = 310 310 310

Pcoupl = c-rescale ; parrinello-rahman
Pcoupltype = isotropic ; semiisotropic
tau_p = 6.0 ; 12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 ; 3e-4
ref_p = 1.0 ; 1.0

gen_vel = yes
gen_temp = 310
gen_seed = -1

constraints = none
constraint_algorithm = Lincs
continuation = no
lincs_order = 4
lincs_warnangle = 30

pbc = xy
nwall = 2
wall-atomtype = W W
wall-density = 10 10
wall-r-linpot = -1.1

NOTE: W is CGMartini water bead.

2

Update:

Even without pressure coupling, after some time GROMACS throws an error saying that an atom is beyond the wall. It happened after around 2 million steps.

3

Update:

I found that decreasing the dt from 0.02 to 0.01 increases the stability. “An atom is beyond the wall” still occurs, but significantly less frequently. Where before it would occur around 2 million steps, with dt = 0.01 I found it to have a ~50% chance to occur around 30,000,000 steps.

4

Maybe the Martini water potential is too soft. Setting wall-r-linpot to a small, positive value might help. You could try 0.1.