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

[Patrick Fuchs]

Hi,
I answer a long time after the last discussion of this thread, but I did 
some tests in between and discussed with some people. I would like to 
follow-up on this issue and sum-up what I found.

After digging into the mailing list, I finally found that the same type 
of issue was discussed in 2011 between Antonio Baptista/Miguel 
Machuqueiro and gromacs devs. There a typical GROMOS setup (twin-range 
cutoff 0.8/1.4 + RF) was giving LINCS crashes (protein simulations in 
water) in 4.5 whereas it was stable in 4.0. It looks like that what we 
see with the Poger lipids has the same origin (the new reversible 
algorithm introduced in 4.5). We don't have crashes because there are no 
explicit H in the lipids. So to sum up the situation, if we want a 
stable simulation with 4.5 and the GROMOS typical simulation parameters, 
that give similar results to 4.0 there seem to be 2 possibilies (maybe 
there are others, but at least these 2 have been tested):
1) reduce nstlist to 1 (2 seems to work according to Miguel)
2) increase rlist up to rlist=rcoulomb (I tried Samuli's suggestion 
rlist=rcoulomb and recover a correct area on the DPPC system, see 
http://redmine.gromacs.org/issues/1400 for the plot).
Both solutions increase a lot the computational cost. For example 2) 
gives 1.5 ns/h on my eight-core computer, whereas in 4.0 I had ~ 1.0 ns/h.

In the discussion of 2011 Berk proposed some other possible solutions to 
speed up with 4.6 
(https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2011-December/005506.html), 
but these were not tested to my knowledge.

My concern is that a twin-range/RF setup in version <= 4.0 gives a given 
result, but does'nt give the same result in >= 4.5 and possibly crashes. 
So I shall say that it would be nice to warn the users that they will 
encounter problems when using twin-range/RF with versions >= 4.5. Maybe 
with a grompp warning? Section 3.4.7 in the manual may not be the first 
place a user will look at.
Anyhow, if for some reason someone wants to use the twin-range cut-off 
as it is done in the GROMOS software (with all versions of GROMOS force 
fields), one has to stick to 4.0.* or lower.

Best,

Patrick

Le 18/12/2013 21:26, Ollila Samuli a écrit :
> Hi,
>
>>> On Dec 18, 2013, at 8:50 AM, Sabine Reisser <sabine.reisser at kit.edu> wrote:
>>> Up to Gromacs 4.0.7, the long rance forces, which are calculated at every nstlist'th step, were added constantly
>>> to the short range force in the subsequent steps. Starting from Gromacs 4.5., the long range force is added
>>> nstlist times to the short range force on the nstlist'th step, while there is nothing added to the short range force
>>> on the subsequent steps. This creates some kind of impulse force which we believe is responsible for the
>>> different behaviour of the membrane, although we haven't figured out the exact mechanism yet.
>
>> That's interesting. Now I see that this change is documented in the manual (Section 3.4.7). Would that imply that >this effect does not depend on the model for long-range electrostatics (Reaction Field or PME), but it comes from >using twin range cutoffs?
>
> I also noticed the same section in the manual today. I think that what possibly happens is that with the reaction field the forces between rlist and rcoulumb are not used in the steps between neighbour list update. However, in these steps the reaction field force calculation assumes the dielectric envinronment beoynd rcoulomb, not beoynd rlist even though the forces only up to rlist are used. I am not sure if this would have any practical relevance, though? I am also not sure how this is handled in the pressure calculation?
>
> I think that the reason why Patrick did not get the difference in the PME results between versions might be that he probably used rcoulomb=rlist in those simulations. If I try to set rcoulomb>rlist with PME, I get error:
>
>    With coulombtype = PME, rcoulomb must be equal to rlist
>    If you want optimal energy conservation or exact integration use PME-Switch
>
> If it is true what I am saying, then the area per molecules should be the same between versions also using RF and setting rlist=rcoulomb.
>
> BR,
> Samuli Ollila
>
>
>
>>> I expect you to get the 'right' area per lipid if you use nstlist 1 as already suggested by Mirco, because in that
>>> case the multistep algorithms give the same answer. This is also our experience with G54A7.
>
>> I started a run yesterday with nstlist=1, and so far (about 12 ns) the area per lipid indeed stays at the "correct" >value of around 0.63 nm^2.
>
> Thanks,
>
>    Lutz
>
> --
> 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.
>

-- 
_______________________________________________________________
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