[gmx-users] mdp-settings for charmm36 and lipid apl values

[Muthukumaran Rajagopalan]

Dear Gromacs Users, 


I have recently started performing membrane simulations in gromacs 5.01 version, I have started running simulations with pure POPE bilayer constructed from charmm-gui membrane builder with XY dimension 70 angs (84 lipids in each bilayer). The bilayer is solvated with 0.15 mol of KCL ( as Na ions have problem with lipid layers). The lipid layer is solvated with hydration number of 36 (lipid:water ratio).


I went through a series of literature and gromacs mailing list regarding decrease in area per lipid and tried different combinations but the area per lipid decreases by 3-4 angs from the experimental values (~58 for POPE). 


I have referred Pluhackova, K et al J. Phys. Chem. B 2016, 120, 3888−3903
Lee et al J. Chem. Theory Comput. 2016, 12, 405−413 

Klauda Jeffery B et al J. Phys. Chem. B, Vol. 114, No. 23, 2010


and previous gromacs mailing list


https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2012-July/073351.html
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2016-June/106385.html



I have used charmm36 force (default given by charmm-gui). I assume charmm-gui provides the charmm specific TIP3P water model (TIPs3P). The simulations start with a minimization with steep and CG methods for (emtol = 100.0), followed by a NVT simulation with restraints on P atoms of lipids ( 10000 KJ/mol) for 200 ps and gradual reduction of restraints 12 steps of NPT each for 500 ps. This is followed by a unrestrained run for 100 ns. I have pasted the md.mdp file options below

integrator = md
dt = 0.002
nsteps = 2500000
nstlog = 5000
nstxout = 0000
nstvout = 0000
nstcalcenergy = 100
nstenergy = 5000
nstxout-compressed = 5000
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
rcoulomb = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 0.8 (as per Pluhackova, K et al J. Phys. Chem. B 2016, 120, 3888−3903 )
rvdw = 1.2

; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
coulomb-modifier = Potential-shift-Verlet (as per Pluhackova, K et al J. Phys. Chem. B 2016, 120, 3888−3903 )
rcoulomb-switch = 0
pme_order = 4 ; cubic interpolation
fourierspacing = 0.12 ; grid spacing for FFT
;

tcoupl = V-rescale
tc_grps = POPE SOL_ION
tau_t = 1.0 1.0
ref_t = 310.00 310.00
;
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 5.0 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
;
; Bond parameters
continuation = yes 
constraint_algorithm = lincs 
constraints = all-bonds 
lincs_iter = 1 
lincs_order = 4 
;
;Periodic boundary conditions
pbc = xyz ; 3-D PBC

; Dispersion correction
DispCorr = no ; account for cut-off vdW scheme

; Velocity generation
gen_vel = no ; Velocity generation is off
 
;
nstcomm = 100
comm_mode = linear
comm_grps = System
;
refcoord_scaling = no

Any suggestions regarding the protocol followed..
 


with regards

Muthukumaran Rajagopalan