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

[Thomas Piggot]

Justin is right about the parameters being sub-optimal, for example see:

de Vries et al. http://pubs.acs.org/doi/abs/10.1021/jp0366926

which I think was the first mention of issues with the 'Berger' PE 
parameters.

There are plenty of force fields out there for PE membranes (both 
all-atom and united-atom) which do a much better job. If you really 
really want to use Berger you can use the fudge as mentioned in this 
paper to include some repulse LJ interactions on the headgroup 
hydrogens. From some tests I did ages ago, I think you need a larger 
repulsive force than is mentioned in the paper (which uses DOPE) to get 
POPE to behave sensibly (IIRC).

Cheers

Tom

On 11/04/15 21:09, Justin Lemkul wrote:
>
>
> On 4/11/15 8:25 AM, niexuechuan wrote:
>> 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?
>>
>
> This could very well be due to suboptimal parameters.  Has anyone ever 
> verified that these POPE parameters produce APL and order parameters 
> in agreement with experimental values over such long time frames?
>
> Lipid force fields are very sensitive and not all parameter sets are 
> created equal.
>
> -Justin
>
>> 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
>>
>

-- 
Dr Thomas Piggot
University of Southampton, UK.