[gmx-users] POPE lipid from Tieleman's web site, the box of xy dimension will shrink during the simulation

[niexuechuan]

Dear all,
   I am new to lipid simulation. I  have downloaded the pope.pdb, pope.itp, lipid.itp from Tieleman's web site (http://wcm.ucalgary.ca/tieleman/downloads). And I have made a directory named gromos53a6_lipid.ff followed Justin's tutorial (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/membrane_protein/index.html). Then I use pope.pdp to run a pure 340 POPE lipid MD at 310K. But I find that the box of xy dimension will shrink from 9.6nm*9.5nm to 9.0nm*8.9nm after 20 ns, then it fluctuates at that value till 200ns. The  box of z dimension will increase a little to hold the volume constant.  This will lead to the area per lipid decreased from 0.54 nm^2 to 0.47nm^2.   At the same time, I find the order parameter of lipid tail (Scd) is too big (above 0.3). Does any one know  where the problem is? 

    When I use dppc128.pdb from Tieleman's web site  to run a pure DPPC lipid MDat 323K, the box will not shrink. And the area per lipid of DPPC is in accordance with literature, about 0.62nm^2. 


  The water model is SPC. My mdp file is also from Justin's tutorial. Followed is my mdp file of POPE:


title= KALP15-DPPC Production MD 
; Run parameters
integrator= md; leap-frog integrator
nsteps= 100000000; 2 * 500000 = 1000 ps (1 ns)
dt    = 0.002; 2 fs
; Output control
nstxout= 0; save coordinates every 2 ps
nstvout= 0; save velocities every 2 ps
nstxtcout= 1000; xtc compressed trajectory output every 2 ps
nstenergy= 1000; save energies every 2 ps
nstlog= 1000; update log file every 2 ps
; Bond parameters
continuation= yes    ; Restarting after NPT 
constraint_algorithm = lincs; holonomic constraints 
constraints= all-bonds        ; all bonds (even heavy atom-H bonds) constrained
lincs_iter= 1            ; accuracy of LINCS
lincs_order= 4            ; also related to accuracy
; Neighborsearching
ns_type= grid; search neighboring grid cels
nstlist= 5    ; 10 fs
rlist= 1.2; short-range neighborlist cutoff (in nm)
rcoulomb= 1.2; short-range electrostatic cutoff (in nm)
rvdw= 1.2; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype= PME; Particle Mesh Ewald for long-range electrostatics
pme_order= 4    ; cubic interpolation
fourierspacing= 0.16; grid spacing for FFT
; Temperature coupling is on
tcoupl= Nose-Hoover    ; More accurate thermostat
tc-grps= POPE SOL;SOL_CL; three coupling groups - more accurate
tau_t= 0.50.5;0.5        ; time constant, in ps
ref_t= 310 310;323        ; 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= 2.0        ; time constant, in ps
ref_p= 1.01.0        ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-54.5e-5; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc    = xyz; 3-D PBC
; Dispersion correction
DispCorr= EnerPres; account for cut-off vdW scheme
; Velocity generation
gen_vel= no; Velocity generation is off
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm         = 1
comm-mode       = Linear
comm-grps       = POPE SOL 
; Scale COM of reference coordinates (when postion restraint)
;refcoord_scaling = com








And my topol file:
#include "gromos53a6_lipid.ff/forcefield.itp" 
#include "pope.itp" 


 ; Include water topology
 #include "gromos53a6_lipid.ff/spc.itp" 


 ; System specifications 
[ system ] 
340-Lipid POPE and CNP
 [ molecules ] 
; molecule name nr. 
POPE 340 
SOL   6729
Thanks in advance,

Nie Xuechuan
 Shanghai Institute of Applied Physics, Chinese Academy of Sciences