Flat-bottom potential renders an unstable system

User discussions
, ,
1

GROMACS version: 2022
GROMACS modification: No

I’m trying to equilibrate a membrane generated with CHARMM-GUI for the Martini forcefield. My membrane contains 512 lipids per monolayer and I’m using 28000 water beads (equivalent to 110 waters per lipid), this membrane as pre-equilibrated using CHARMM-GUI’s protocol. I want to make a pore in the membrane using the Flat-Bottom potential to keep the lipid tails out of a region on the z-axis. I’m following the CHARMM-GUI’s protocol for equilibrating vesicles, but when I apply the potential, the system symply crashes, most of the time mdrun just freezes and refuses to run. It won’t even give any error except sometimes a segmentatium fault

I have the following lines in my itp file under the lipids definition:

;keep lipid tails out of a cylinder along X Y Z axis to maintain the waterpore:
#ifdef  VESICLE_LIPIDTAIL_R
    #ifndef VESICLE_LIPIDTAIL_FC
        #define VESICLE_LIPIDTAIL_FC 5
    #endif
    [ position_restraints ]
      5     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
      6     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
      7     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
      8     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
      9     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
     10     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
     11     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
     12     2    2   -VESICLE_LIPIDTAIL_R  VESICLE_LIPIDTAIL_FC
#endif

and my mdp file is as follows:

define                   = -DVESICLE_LIPIDTAIL_R=2.00 
integrator               = md
tinit                    = 0.0
dt                       = 0.002
nsteps                   = 1000000

nstlog                   = 1000
nstenergy                = 1000
nstxout-compressed       = 1000
compressed-x-precision   = 100

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

epsilon_r                = 15
coulombtype              = reaction-field
rcoulomb                 = 1.1
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 1.1

tcoupl                   = v-rescale
tc-grps                  = membrane solute
tau_t                    = 1.0  1.0
ref_t                    = 303.15 303.15

; Pressure coupling:
Pcoupl                   = berendsen
Pcoupltype               = semi-isotropic
tau_p                    = 5.0
compressibility          = 4.5e-5   4.5e-5
ref_p                    = 1.0  1.0

; GENERATE VELOCITIES FOR STARTUP RUN:
gen_vel                  = yes
gen_temp                 = 303.15
gen_seed                 = 8415458562

refcoord_scaling         = all

If I decreade the -DVESICLE_LIPIDTAIL_R value to 0.50 or lower the simulation runs, but it won’t produce the desired effect. Any value of #define VESICLE_LIPIDTAIL_FC greater than 5 also makes the simulation crash

How can I avoid the system crashing?? I thought perhaps my system is too small, but a bilayer with 512 lipids per monolayer is already pretty big. So I’m clueless on how to make this work

Thanks