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

[Justin Lemkul]

On Mon, Dec 16, 2013 at 3:10 PM, Lutz Maibaum <lutz.maibaum at gmail.com>wrote:

> 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.
>
>
I suspect you're running into the same issue that has been reported here:
http://redmine.gromacs.org/issues/1400.

Can you check to see if Reaction-Field-nec fixes the issue?  We're still
waiting on feedback from the Redmine issue.  The RF methods changed
somewhere along the way, and we need to make sure that the appropriate
algorithms are being tested for an apples-to-apples comparison.

-Justin

-- 

==========================================

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441


==========================================