[gmx-users] Using CHARMM 36 for DPPC simulation

[Amit Choubey]

5Hello all,

I am trying to use CHARMM 36 for DPPC membrane simulation. I did the
following so far:

1. Download pdb file containing 128 DPPC molecules from
http://www.charmm-gui.org/?doc=archive&lib=lipid_pure

2. I separated one lipid molecule from the obtained pdb file and used
pdb2gmx -f 1dppc.gro -nochargegrp

The top file obtained was converted to itp file by commenting few things
out.

3. I looked at a file named dppc_n128.inp in the downloaded dppc bilayer
from CHARMM GUI for the box length. I used editconf to convert the pdb to
gro and manually inserted the box size in the gro file.

4. Created a system.top file which looks like
#include "charmm36.ff/forcefield.itp"
#include "lipid.itp"
#include "charmm36.ff/tips3p.itp"

[ system ]
chol

[ molecules ]
DPPC 128
SOL  3659

5. I used following mdp settings

integrator = md ; leap-frog integrator
nsteps = 20000000 ; 40 ns
dt = 0.002 ; 2 fs

; Output control
nstxout = 100000 ; save coordinates every 200 ps
nstvout = 100000 ; save velocities every 200 ps
nstenergy = 10000 ; save energies every 2 ps
nstlog = 10000 ; update log file every 2 ps

; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 10 ; 10 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
rlistlong = 1.4
vdwtype = switch
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
rvdw_switch = 0.8

; Electrostatics
rcoulomb        = 1.0
rcoulomb_switch = 0.0
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 6 ; cubic interpolation
fourierspacing = 0.15 ; grid spacing for FFT

; Temperature coupling is on
tcoupl = Nose-Hoover ; More accurate thermostat
tc-grps = DPPC SOL  ; two coupling groups - more accurate
tau_t = 0.2 0.2 ; time constant, in ps
ref_t = 323.15 323.15 ; reference temperature, one for each group, in K

; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors,
independent z
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1

; Periodic boundary conditions
pbc = xyz ; 3-D PBC

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


; Velocity generation
gen_vel             =  yes
gen_temp            =  200.0
gen_seed            =  173529

; COM motion removal
nstcomm         = 1
comm-mode       = Linear
comm-grps       = DPPC SOL

; Energy monitoring
energygrps          =  DPPC SOL

; Bond parameters
constraint_algorithm = lincs    ; holonomic constraints
constraints     = hbonds     ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy

The resulting simulation gave me an area per lipid =0.591 nm^2. This is not
quite right. Can somebody suggest what else could i do it get close to the
0.63 nm^2 value.

Also when i looked at trajectory file it seemed that the box length was not
set up quite right because the downloaded pdb file looked ok with the
membrane in the middle but at the next frame the membrane was at a totally
different location. The leaflets were separated and the middle of the
membrane was at the edge of the box.

Also when i used pdb2gmx -nochargegrp on the downloaded file it created a
gro file whose system size was
   8.20170   8.64120   6.54740 in contrast to    6.40000   6.40000
6.38452 .

Amit
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111109/ae4089f7/attachment.html>