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

[Patrick Fuchs]

Hi,
to follow up on this, the simulation with "Reaction-Field-nec" under 
4.5.3 has completed. The final area is 0.60 nm^2 (see 
http://redmine.gromacs.org/issues/1400 for the plot). Lutz, you 
mentionned you tried a simulation with nstlist=1, could you give 
feedback on redmine?
Best,

Patrick

Le 16/12/2013 23:14, Patrick Fuchs a écrit :
> 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