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

[Patrick Fuchs]

Hi,
following Justin's suggestion on redmine (issue #1400), a simulation 
with "Reaction-Field-nec" using 4.5.3 is running. It hasn't completed 
yet but the mean area per lipid around 50 ns is 0.59 < mean area < 0.60 
nm^2. I will update on redmine once the simulation is completed 
(probably in a couple of days).
Best,

Patrick

Le 16/12/2013 22:30, Justin Lemkul a écrit :
> 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
>

-- 
_______________________________________________________________
Patrick FUCHS
Dynamique des membranes et trafic intracellulaire
Institut Jacques Monod, CNRS UMR 7592, Université Paris Diderot
Bâtiment Buffon, 15 rue Hélène Brion, 75013 Paris
Tel :  +33 (0)1 57 27 80 05 - Fax : +33 (0)1 57 27 81 35
E-mail address: patrick.fuchs at univ-paris-diderot.fr
Web Site: http://www.dsimb.inserm.fr/~fuchs