[gmx-users] mdp-settings for charmm36 and lipid apl values
Christopher Neale
chris.neale at mail.utoronto.ca
Thu May 9 01:06:09 CEST 2013
It is a really bad idea to use standard tip3p with charmm36 lipids (see the Piggot paper that you referenced and also Sapay, N. et al. 2010 J. Comp. Chem. 32, 1400-1410 + probably others).
dt 0.001 with nstlist 5 seems like overkill on the nstlist update frequency (not a problem though).
Here's how I do charmm36 lipid simulations:
constraints = all-bonds
lincs-iter = 1
lincs-order = 6
constraint_algorithm = lincs
integrator = sd
dt = 0.002
tinit = 0
nsteps = 5000000
nstcomm = 1
nstxout = 5000000
nstvout = 5000000
nstfout = 5000000
nstxtcout = 50000
nstenergy = 50000
nstlist = 10
nstlog=0 ; reduce log file size
ns_type = grid
vdwtype = switch
rlist = 1.2
rlistlong = 1.3
rvdw = 1.2
rvdw-switch = 0.8
rcoulomb = 1.2
coulombtype = PME
ewald-rtol = 1e-5
optimize_fft = yes
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
tc_grps = System
tau_t = 1.0
ld_seed = -1
ref_t = 310
gen_temp = 310
gen_vel = yes
unconstrained_start = no
gen_seed = -1
Pcoupl = berendsen
pcoupltype = semiisotropic
tau_p = 4 4
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
For a pure POPC bilayer with Charmm36 lipids and tips3p water, I get ~ 0.64 nm^2/lipid (not using grid-mat, just looking at box dimensions).
If I use regular tip3p instead, the APL decreases a lot and eventually forms a gel. In simple terms of APL and phase, you can get the same results as tips3p (0.64 nm^2/lipid) if you use tip4p (and spce is not too bad, but is still not as close as tip4p). Note tip4p will run faster than tips3p by a large margin. Nevertheless, I use tip3sp in all of my work with charmm36.
Also, note that ions can shrink your APL, particularly Na+ (and divalent cations are even worse). My simulations use 50 mM KCl, which basically doesn't affect the average APL (although 50 mM NaCl does noticeably reduce the APL).
Finally, I am not convinced that a per-molecule area per lipid is a useful quantity to compare to experimental areas per lipid. I haven't looked at Grid-Mat myself, but there must be a lot of assumptions underlying any analysis that tried to assign an "area" to a single 3D lipid. If I were you, I'd be looking for expt. results of APL for POPC vs. POPE vs. 70%POPC/30%POPE. Also, are you entirely sure that you didn't mix up the POPC vs. POPE values? It looks to me like a simple labelling error enough that it warrants a second look.
Chris.
-- original message --
I've been experimenting with simulations of mixed bilayers (512 lipids in
total, 70% POPC, 30% POPE) using the charmm36 parameter set in gromacs, and
have a couple of questions. I know this has been discussed before, but I'd
appreciate some input nonetheless :-)
The relevant sections of my mdp-file are pasted below:
; Start time and timestep in ps
tinit = 0
dt = 0.001
nsteps = 100000000
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 5
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic_molecules = no
; nblist cut-off
rlist = 1.2
; long-range cut-off for switched potentials
rlistlong = 1.4
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 1.2
; Relative dielectric constant for the medium and the reaction field
epsilon_r = 1
epsilon_rf = 1
; Method for doing Van der Waals
vdw-type = switch
; cut-off lengths
rvdw-switch = 0.8
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = No
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = V-rescale
nsttcouple = -1
nh-chain-length = 10
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.1
ref_t = 300
; Pressure coupling
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
nstpcouple = -1
This is as far as I can tell from earlier discussions on the list, and also
from reading the Piggott et al. paper in JCTC, the correct settings for
charmm36.
After a simulation of ~50 ns, I use GridMatMD to calculate the area per
headgroup of POPC and POPE, respectively, and get what I think are not 100%
acceptable results (but maybe they are)
For POPC, I get 59,7 A^2, and for POPE, I get 63,1 A^2.
The value for POPE would have been fine I suppose if it hadn't been for the
fact that the APL for POPC is smaller. Should it not be larger than POPE?
I notice in the Piggott-paper that they in the supplement for some
simulations of POPC also get APL's of around 59-60 (without POPE of
course), and that the results depend to some extent also on the usage of
TIP3P vs TIP3SP water models. I have been using "normal" TIP3P here.
Could anyone comment on (a) my mdp-file settings, and (b) the resulting
APL, and tell me if I should be worried about anything?
THANKS
More information about the gromacs.org_gmx-users
mailing list