[gmx-users] gromacs (git) with charmm FF lipid freezing problem

[David Caplan]

Hello,

I am trying to simulate a DPPC bilayer, 180 lipids total, ~13000 TIP3
and 50 K, 50 Cl counter-ions. The standard charmm27r DPPC parameters
have very large charge groups, so I went through them and made them
smaller (4 atoms each, on average), following the comments from
grompp. I have tried the C27 parameters, recent parameters by Sonne,
et al. (Reparameterization of All-Atom Dipalmitoylphosphatidylcholine
Lipid Parameters Enables Simulation of Fluid Bilayers at Zero Tension.
Jacob Sonne, et al. Biophysical Journal, 15 June 2007. Vol. 92, Issue
12, pp. 4157-4167) as well as the latest Charmm36 parameters available
from Alex MacKerrell's website. The last two do not require surface
tension so I was using semiisotropic pressure coupling with the same
compressibility for XY and Z (4.5e-5). For charmm27 parameters I set
XY compressibility to 0 (constant area).

I am using all-bond constraints with LINCS and the CHARMM_TIP3P water
model. I have tried langevin dynamics as well as the standard md
integrator. The temperature is coupled separately for DPPC and
non-DPPC at 323.15K.

When I run any simulation for >10 ns, I get what looks like frozen
lipids. They become very straight and aligned with each other,
sometimes intercalated. I'm not sure why this is happening... It may
be that the lipids are using incorrect dihedral and/or angle
parameters (only charges are set explicitly in my DPPC ITP files,
everything else is in the base ffnb.itp)

Any guesses as to why this may be happening? Anything specific I
should be analyzing? I see a steady decrease in LR coulombic potential
which seems to change depending on my fourierspacing... so I wonder if
my PME settings are wrong.

Here is part my MDP:

integrator               = sd
define  = -DCHARMM_TIP3P
nsteps = 500000000
dt = 0.002
bd-fric                  = 0
ld-seed                  = -1
nstxout         = 0
nstvout         = 0
nstfout         = 0
nstlog          = 500
nstenergy       = 500
nstxtcout       = 1000
xtc-precision   = 1000
nstcomm = 5
nstcalcenergy = 5
nstlist     = 10
ns_type     = grid
pbc     = xyz
coulombtype = PME
rcoulomb = 1.3
pme_order = 6
ewald_rtol = 1e-5
vdwtype = Switch
rvdw            = 1.3
rvdw_switch     = 1.0
rlist           = 1.3
rlistlong       = 1.5
fourierspacing = 0.12
tc-grps         = DPPC non-DPPC
tau_t            = 1 1
ref_t             = 323.15 323.15
pcoupl                  = Berendsen
pcoupltype               = semiisotropic
tau_p                    = 1.0
ref_p                    = 1.0 1.0
compressibility          = 4.5e-5 4.5e-5
gen_vel                  = yes
gen_temp                 = 323
gen_seed                 = -1
constraints              = all-bonds
constraint-algorithm     = Lincs
lincs-iter  = 1
lincs-order              = 6
lincs-warnangle          = 60
morse                    = no