Charged groups penetrating into hydrophobic region

GROMACS version: 2022.4
GROMACS modification: No

Hi,

I’m trying to simulate a POPG membrane with 128 lipids per monolayer and a pore containing 12 peptides. But every time I run an energy minimization the phopholipids phosphate groups and even the Na+ conterions start penetrating the membrane hydrophobic core.

I’ve check the system net charge, and it’s equal to zero, Gromacs also didn’t give me any note or warning about the system charge.

Most importantly the charged groups start penetrating in the entire membrane, not on the pore region.
I’ve also tried to run a membrane without the pore with the same files but I get the same problem.
I’ve tried to do the minimization and equilibration with restraints, but it doesn’t work the same problem occurs when I remove the restraints.

I assembled my membrane with the pore with the CHARMM-GUI bilayer builder.
I’m using Slipids for the membrane and Amber ff99sb-ildn-nmr for the peptides .
I downloaded the .itp file for the POPG lipid from the Slipids web site.

The membrane starts as it’s on the left and aftyer the minimization gets to what is on the right. The phosphates are shown in orange and the sodium ions in purple.

my energy minimization .mdp file is as follows:

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save

define          =  -DPOSRES ;Position restraints for the peptides
integrator      = steep         ; Algorithm (steep = steepest descent minimization)
emtol           = 1        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nmi ;  the minimization converges when the max force is smaller than this value (in units of kJ mol –1 nm –1)
emstep          = 0.001          ; Energy step size
nsteps          = 50000          ; Maximum number of (minimization) steps to perform

cutoff-scheme   = Verlet
ns_type         = grid          ; search neighboring grid cells
vdw_modifier    = potential_switch
rcoulomb        = 1.5           ; short-range electrostatic cutoff (in nm)
rvdw            = 1.5           ; short-range van der Waals cutoff (in nm)
rvdw_switch     = 1.4

; Electrostatics
coulombtype     = PME           ; Particle Mesh Ewald for long-range electrostatics
pme_order       = 4             ; cubic interpolation
fourierspacing  = 0.16          ; grid spacing for FFT

All my files can be found on this link: Google Drive

Any help will be really apreciated.

Thanks