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

[Ollila Samuli]

Hi,

Essentially similar results have been reported recently in redmine: http://redmine.gromacs.org/issues/1400

The issue seems to be related to the reaction field since with PME the results were the same between the versions. It might be useful to check if you get the similar pressure differences between versions also with PME? If not, then the pressure calculation with RF in 4.5. and 4.6 might be the issue. In this case also the different components of pressure might be interesting.

Maybe it is better to report this and further findings by updating redmine issue 1400?

BR,
Samuli Ollila


________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Lutz Maibaum [lutz.maibaum at gmail.com]
Sent: Monday, December 16, 2013 10:10 PM
To: gromacs.org_gmx-users at maillist.sys.kth.se
Subject: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and       4.6.5

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)



--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.