<div><br></div>5Hello all,<div><br></div><div>I am trying to use CHARMM 36 for DPPC membrane simulation. I did the following so far:</div><div><br></div><div>1. Download pdb file containing 128 DPPC molecules from <a href="http://www.charmm-gui.org/?doc=archive&lib=lipid_pure">http://www.charmm-gui.org/?doc=archive&lib=lipid_pure</a></div>
<div><br></div><div>2. I separated one lipid molecule from the obtained pdb file and used </div><div>pdb2gmx -f 1dppc.gro -nochargegrp </div><div><br></div><div>The top file obtained was converted to itp file by commenting few things out.</div>
<div><br></div><div>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. </div>
<div><br></div><div>4. Created a system.top file which looks like</div><div><div>#include "charmm36.ff/forcefield.itp"</div><div>#include "lipid.itp"</div><div>#include "charmm36.ff/tips3p.itp"</div>
<div><br></div><div>[ system ]</div><div>chol</div><div><br></div><div>[ molecules ]</div><div>DPPC 128</div><div>SOL 3659</div></div><div><br></div><div>5. I used following mdp settings</div><div><br></div><div><div>integrator<span class="Apple-tab-span" style="white-space:pre"> </span>= md<span class="Apple-tab-span" style="white-space:pre"> </span>; leap-frog integrator</div>
<div>nsteps<span class="Apple-tab-span" style="white-space:pre"> </span>= 20000000<span class="Apple-tab-span" style="white-space:pre"> </span>; 40 ns</div><div>dt<span class="Apple-tab-span" style="white-space:pre"> </span>= 0.002<span class="Apple-tab-span" style="white-space:pre"> </span>; 2 fs</div>
<div><br></div><div>; Output control</div><div>nstxout<span class="Apple-tab-span" style="white-space:pre"> </span>= 100000<span class="Apple-tab-span" style="white-space:pre"> </span>; save coordinates every 200 ps</div>
<div>nstvout<span class="Apple-tab-span" style="white-space:pre"> </span>= 100000<span class="Apple-tab-span" style="white-space:pre"> </span>; save velocities every 200 ps</div><div>nstenergy<span class="Apple-tab-span" style="white-space:pre"> </span>= 10000<span class="Apple-tab-span" style="white-space:pre"> </span>; save energies every 2 ps</div>
<div>nstlog<span class="Apple-tab-span" style="white-space:pre"> </span>= 10000<span class="Apple-tab-span" style="white-space:pre"> </span>; update log file every 2 ps</div><div><br></div><div>; Neighborsearching</div>
<div>ns_type<span class="Apple-tab-span" style="white-space:pre"> </span>= grid<span class="Apple-tab-span" style="white-space:pre"> </span>; search neighboring grid cels</div><div>nstlist<span class="Apple-tab-span" style="white-space:pre"> </span>= 10<span class="Apple-tab-span" style="white-space:pre"> </span>; 10 fs</div>
<div>rlist<span class="Apple-tab-span" style="white-space:pre"> </span>= 1.0<span class="Apple-tab-span" style="white-space:pre"> </span>; short-range neighborlist cutoff (in nm)</div><div>rlistlong <span class="Apple-tab-span" style="white-space:pre"> </span>= 1.4</div>
<div>vdwtype<span class="Apple-tab-span" style="white-space:pre"> </span>= switch</div><div>rvdw<span class="Apple-tab-span" style="white-space:pre"> </span>= 1.2<span class="Apple-tab-span" style="white-space:pre"> </span>; short-range van der Waals cutoff (in nm)</div>
<div>rvdw_switch<span class="Apple-tab-span" style="white-space:pre"> </span>= 0.8</div><div><br></div><div>; Electrostatics</div><div>rcoulomb = 1.0</div><div>rcoulomb_switch<span class="Apple-tab-span" style="white-space:pre"> </span>= 0.0</div>
<div>coulombtype<span class="Apple-tab-span" style="white-space:pre"> </span>= PME<span class="Apple-tab-span" style="white-space:pre"> </span>; Particle Mesh Ewald for long-range electrostatics</div><div>pme_order<span class="Apple-tab-span" style="white-space:pre"> </span>= 6<span class="Apple-tab-span" style="white-space:pre"> </span>; cubic interpolation</div>
<div>fourierspacing<span class="Apple-tab-span" style="white-space:pre"> </span>= 0.15<span class="Apple-tab-span" style="white-space:pre"> </span>; grid spacing for FFT</div><div><br></div><div>; Temperature coupling is on</div>
<div>tcoupl<span class="Apple-tab-span" style="white-space:pre"> </span>= Nose-Hoover<span class="Apple-tab-span" style="white-space:pre"> </span>; More accurate thermostat</div><div>tc-grps<span class="Apple-tab-span" style="white-space:pre"> </span>= DPPC<span class="Apple-tab-span" style="white-space:pre"> </span>SOL ; two coupling groups - more accurate</div>
<div>tau_t<span class="Apple-tab-span" style="white-space:pre"> </span>= 0.2<span class="Apple-tab-span" style="white-space:pre"> </span>0.2 ; time constant, in ps</div><div>ref_t<span class="Apple-tab-span" style="white-space:pre"> </span>= 323.15 <span class="Apple-tab-span" style="white-space:pre"> </span>323.15 ; reference temperature, one for each group, in K</div>
<div><br></div><div>; Pressure coupling is on</div><div>pcoupl<span class="Apple-tab-span" style="white-space:pre"> </span>= Parrinello-Rahman<span class="Apple-tab-span" style="white-space:pre"> </span>; Pressure coupling on in NPT</div>
<div>pcoupltype<span class="Apple-tab-span" style="white-space:pre"> </span>= semiisotropic<span class="Apple-tab-span" style="white-space:pre"> </span>; uniform scaling of x-y box vectors, independent z</div><div>tau_p<span class="Apple-tab-span" style="white-space:pre"> </span>= 5.0<span class="Apple-tab-span" style="white-space:pre"> </span>; time constant, in ps</div>
<div>ref_p<span class="Apple-tab-span" style="white-space:pre"> </span>= 1.0<span class="Apple-tab-span" style="white-space:pre"> </span>1.0<span class="Apple-tab-span" style="white-space:pre"> </span>; reference pressure, x-y, z (in bar)</div>
<div>compressibility = 4.5e-5 4.5e-5<span class="Apple-tab-span" style="white-space:pre"> </span>; isothermal compressibility, bar^-1</div><div><br></div><div>; Periodic boundary conditions</div><div>pbc<span class="Apple-tab-span" style="white-space:pre"> </span>= xyz<span class="Apple-tab-span" style="white-space:pre"> </span>; 3-D PBC</div>
<div><br></div><div>; Dispersion correction</div><div>DispCorr<span class="Apple-tab-span" style="white-space:pre"> </span>= No<span class="Apple-tab-span" style="white-space:pre"> </span>; account for cut-off vdW scheme</div>
<div><br></div><div><br></div><div>; Velocity generation</div><div>gen_vel = yes</div><div>gen_temp = 200.0</div><div>gen_seed = 173529</div><div><br></div><div>; COM motion removal</div>
<div>nstcomm = 1</div><div>comm-mode = Linear</div><div>comm-grps = DPPC SOL</div><div><br></div><div>; Energy monitoring</div><div>energygrps = DPPC SOL</div><div><br></div><div>; Bond parameters</div>
<div>constraint_algorithm = lincs ; holonomic constraints</div><div>constraints = hbonds ; all bonds (even heavy atom-H bonds) constrained</div><div>lincs_iter = 1 ; accuracy of LINCS</div><div>
lincs_order = 4 ; also related to accuracy</div></div><div><br></div><div>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. </div>
<div><br></div><div>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.</div>
<div><br></div><div>Also when i used pdb2gmx -nochargegrp on the downloaded file it created a gro file whose system size was </div><div> 8.20170 8.64120 6.54740 in contrast to 6.40000 6.40000 6.38452 . </div>
<div><br></div><div>Amit</div>