[gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5

[Lutz Maibaum]

I am running MD simulations of DPPC bilayers with the Gromos force field, and I am seeing some differences between using Gromacs 4.0.7 and 4.5.5/4.6.5 that I do not understand. It would be great if someone had any insight into what's going on here.

When I use the Gromos 53A6-L force field, which is the default 53A6 force field combined with improved lipid parameters that can be downloaded from http://compbio.chemistry.uq.edu.au/~david/research/lipids.htm, I obtain an average area per lipid of 0.627 nm^2, in good agreement with both the paper that describes these new parameters (0.629 nm^2, obtained with Gromacs 3.2.1, Ref. [1]) and another follow-up study (0.631 nm^2 and 0.623 nm^2, Gromacs 4.0.7, Ref. [2]).

Now, if I use the same force field and mdp file, and the same initial configuration (which is a pre-equilibrated DPPC bilayer from http://compbio.biosci.uq.edu.au/atb/system_download.py?boxid=32 and randomly generated velocities), but use Gromacs 4.6.5 instread, I get a lower value of about 0.59 nm^2.

I also tried the Gromos 54A7 force field, which is included with Gromacs 4.6.5 and that should be identical to 53A6-L (plus it has some other improvements over 53A6 that shouldn't be relevant here), I also get the lower area per lipid of ~0.59 nm^2.

If have attached a plot of the area per lipid for these three simulations, each more than 100ns long. This looks to me like 4.6.5 give systematically lower area per lipid than 4.0.7. Running additional simulations with Gromacs 4.5.5 suggest that that also results in the lower area per lipid. Does anyone know why this might be?  I have uploaded the relevant files in case that is helpful:

http://faculty.washington.edu/maibaum/dppc_comparison/

To see if there are any differences between the energies that 4.0.7 and 4.6.5 compute, I picked a configuration, and used "mdrun -rerun" with the three different gromacs / force field combinations. Here is what I get for the "sample.gro" configuration (included in the link above):

Gromacs 4.0.7 + Gromos53A6-L:

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    1.76906e+02    1.29521e+04    9.57922e+03    4.97638e+02   -1.38126e+03
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coulomb (LR)
    1.50791e+04    1.23533e+04   -7.49664e+03   -3.78489e+05   -8.97363e+02
       RF excl.      Potential    Kinetic En.   Total Energy    Temperature
   -5.25414e+04   -3.90168e+05    6.47185e+04   -3.25449e+05    2.87013e+02
 Pressure (bar)
    2.00740e+02


Gromacs 4.6.5 + Gromos53A6-L:

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    1.76905e+02    1.29522e+04    9.57922e+03    4.97637e+02   -1.38126e+03
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coulomb (LR)
    1.50791e+04    1.23533e+04   -7.49659e+03   -3.78489e+05   -8.97344e+02
       RF excl.      Potential    Kinetic En.   Total Energy    Temperature
   -5.25414e+04   -3.90168e+05    6.47458e+04   -3.25422e+05    2.87134e+02
 Pressure (bar)
    2.15012e+02


Gromacs 4.6.5 + Gromos54A7:

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    1.76905e+02    1.30332e+04    9.57922e+03    4.97637e+02   -1.38126e+03
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coulomb (LR)
    1.50791e+04    1.25510e+04   -7.32610e+03   -3.78489e+05   -8.97344e+02
       RF excl.      Potential    Kinetic En.   Total Energy    Temperature
   -5.25414e+04   -3.89718e+05    6.47631e+04   -3.24955e+05    2.87211e+02
 Pressure (bar)
    2.89622e+02


I don't see any significant difference between what 4.0.7 and what 4.6.5 compute. I don't know if the somewhat higher pressure with the 4.6.5/54A7 combination is meaningful.

If anyone can shed any light on this, or has ideas for how to debug this further, I'd be most grateful.

Best,

  Lutz


[1] A new force field for simulating phosphatidylcholine bilayers. D. Poger and W. F. van Gunsteren and A. E. Mark. J. Comp. Chem.  31  1117--1125  (2010)

[2] Molecular Dynamics Simulations of Phosphatidylcholine Membranes: A Comparative Force Field Study. T. J. Piggot and A. Pineiro and S. Khalid. J. Chem. Theo. Comp.  8  4593--4609  (2012)