[gmx-users] Gromos54a7 electrostatics interactions ?

[Thomas Piggot]

On 22/09/17 16:59, Sim gmx wrote:
> First of all, thanks for your replies!
>
> 1) Indeed, both epsilon-r and epsilon-rf are such constants. However, in
> the paper I guess that they refer to the relative dielectric constant of
> the reaction field, i.e. epsilon-rf. Or do I miss something ?
Yes, they are referring to epsilon-rf.
> 2) I will for sure contact them to collect all these info if I go further
> with this forcefield!
I've tested things pretty extensively with this force field so I don't 
think you need to contact them. If you post your mdp it will make things 
much easier as we can see exactly what you're are doing and can let you 
know if there is anything that doesn't look correct.
>
> 3) I don't think I understand everything there. What I think I get is:
> using a verlet scheme with parameters that were not initially designed for
> this scheme is risky and should be carefully assessed.
This is what I was saying I have already assessed and, for DPPC at 
least, the membrane does not behave well with Verlet and reaction-field. 
If you use PME, Verlet and group give the same results (but the membrane 
is more ordered with PME compared to RF, irrespective of the cutoff scheme).
>   Plus, parameters
> used by Poger are very likely to be incompatible with verlet scheme.
>  From this, I conclude it's not even worth giving it a try. So my options
> are:
> - use a group scheme, that will be slower
Well, this depends upon your computer. On one HPC I use that is CPU 
only, I get nearly identical performance with Verlet and a straight 1.4 
nm cutoff with RF (i.e the settings that don't work well) and the 
fastest group twin-range setup that gives appropriate membrane 
properties. If you are wanting to use a GPU this will obviously be 
different.
> - use a different force field, but which one since I do have the topologies
> of my small molecules for gromos 53a6/54a7 ? Would Berger Lipids be a good
> option here ?
This really depends upon what you are wanting to do. There are other 
GROMOS compatible lipid force fields out there (including Berger) but, 
as with pretty much everything, there are limitations/considerations 
with them all. In addition, you say you already have small molecule 
topologies. Depending upon what these structures are, it may be fairly 
trivial to make ones for other force fields too. This would open up lots 
of other possible lipid force fields too. Again depending upon what you 
are wanting to look at, these could be better/more appropriate for your 
study. Personally, I'd suggest looking at more than once force field, 
that way you can be more confident in any results you obtain. This, 
obviously, takes more time and effort though.

You should probably do a bit of reading about the 
advantages/disadvantages of different (lipid) force fields, that way you 
could make your own mind up rather than just relying on someones opinion 
on a mailing list! Although if you do want my opinion, I'm happy to 
share it if you can give more details of what you are trying to do.
>
> About Tom's reply:
>
> 1. That's good to know, especially as it does not appear anywhere in the
> original papers.
The original simulations in the papers were performed with a really old 
version of GROMACS (3.2.1 IIRC) and so this was not an issue. The paper 
(by the same authors) I linked to in my previous message provides all 
the details about this.
> 2. Is using a setup with all cutoffs = 1.4 more "adequate" than the
> twin-range setup then ? Or is it just a "good to know" stuff ? I would
> naturally use the fastest correct setup which seems to be the twin-range
> one (with group-scheme), am I wrong ?
With the group cutoff scheme, simulations using both the twin-range 
0.8/1.4 and the straight 1.4 cutoffs result in the same properties of a 
DPPC membrane (again see the paper I linked to; I have also tested this 
too). So yes, you should use the faster setup.
>
> Again, thanks to both of you.
>
> 2017-09-22 14:29 GMT+02:00 Piggot T. <T.Piggot at soton.ac.uk>:
>
>> Hi,
>>
>> In addition to Mark's comments, there are a few other points to be aware
>> of:
>>
>> 1. If you use the twin-range scheme (i.e. 0.8/1.4) with nstlist 5, you
>> need to set nstcalclr to 1 to get results matching those as reported by
>> Poger and co-workers (see http://pubs.acs.org/doi/abs/
>> 10.1021/acs.jctc.7b00178?src=recsys). This is the fastest setup that will
>> give you the 'correct' membrane properties with the group scheme.
>>
>> 2. The combination of reaction field with Verlet and all cutoffs set to
>> 1.4 does not give the correct membrane properties for DPPC (these are my
>> own results). Using the group scheme with this setup is fine but slower
>> than the twin-range setup mentioned above. You should therefore use the
>> group setup if using reaction-field with this force field. If you wish to
>> use PME the membrane will be more ordered. However, Verlet and group
>> schemes give the same results for PME (again my results),
>>
>> Regarding the second point, can anyone shed some more light onto why using
>> the Verlet scheme would result in substantially different membrane
>> properties to the group scheme with a single 1.4 nm cutoff and
>> reaction-field (PME gives the same membrane properties using these two
>> setups)?
>>
>> Cheers
>>
>> Tom
>> ________________________________________
>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [
>> gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Mark
>> Abraham [mark.j.abraham at gmail.com]
>> Sent: 22 September 2017 12:41
>> To: gmx-users at gromacs.org; gromacs.org_gmx-users at maillist.sys.kth.se
>> Subject: Re: [gmx-users] Gromos54a7 electrostatics interactions ?
>>
>> Hi,
>>
>> On Fri, Sep 22, 2017 at 10:21 AM Sim gmx <simgmx at gmail.com> wrote:
>>
>>> Hi!
>>>
>>> I would like to do simulations of lipids bilayers with gromos54a7
>>> parameters. To do so, I want to use the same mdp parameters as Poger et
>> al.
>>> (the authors).
>>>
>>> In their papers, they write : "Nonbonded interactions were evaluated
>> using
>>> a twin-range cutoff scheme: interactions falling within the 0.8-nm
>>> short-range cutoff were calculated every step, whereas interactions
>> falling
>>> within the 1.4-nm long-range cutoff were updated every 10 fs, together
>> with
>>> the pair list. A reaction-field correction was applied to account for the
>>> truncation of the electrostatic interactions beyond the long-range cutoff
>>> using a relative dielectric permittivity constant of 62, as appropriate
>> for
>>> SPC water"
>>>
>> Note that such twin-range multiple stepping schemes have only ever been
>> supported in GROMACS in the group cutoff scheme, and not at all in the most
>> recent versions. See (older versions of) the reference manual for
>> discussion of known issues with the twin-range scheme. However, parameters
>> derived with this irreversible (and thus unphysical) scheme have usually
>> been shown to be usable outside that scheme, but obviously there's a higher
>> burden of validation bourne by a the user in that case.
>>
>>  From this, I conclude that my mdp file should contain something like (among
>>> other parameters):
>>>
>>> dt=0.002
>>> vdw-type = cut-off
>>> rlist = 0.8
>>> nstlist = 5
>>> rvdw = 1.4
>>> coulombtype = reaction-field
>>> epsilon-rf = 62
>>>
>>> My questions are:
>>>
>>> 1) Does "relative dielectric permittivity constant" indeed mean
>>> "epsilon-rf" here (not epsilon-r) ?
>>>
>> That's not the right question, both epsilon and epsilon-rf are such
>> constants. See
>> http://manual.gromacs.org/documentation/2016.4/user-
>> guide/mdp-options.html#electrostatics
>> for
>> the different parts of the coulomb contributions to which they refer.
>>
>>
>>> 2) What about the coulomb cut-off (rcoulomb) ? I think that in
>> gromos53a6,
>>> rcoulomb is often set to 0.9 (here I guess it would be 0.8 to match the
>>> rlist). However, if the reaction-field correction is applied "to account
>>> for the truncation of the electrostatic interactions beyond *the
>> long-range
>>> cutoff*", I would set rcoulomb to 1.4. As no explicit reference to
>> rcoulomb
>>> is done in the papers, I fear to miss something here.
>>>
>> They refer to non-bonded interactions, that means both Coulomb and VDW.
>> Your interpretation is correct for the group scheme. Their text implies
>> rlist=0.8 and the others 1.4. You should also ask them if they used GROMACS
>> and to supply a sample .mdp file so you can reproduce their reported
>> results.
>>
>> 3) If rcoulomb has indeed to be set to 1.4 as rvdw, it becomes possible to
>>> use a verlet cut-off scheme to benefit from GPU acceleration. However,
>> when
>>> verlet cut-off scheme is used, mdrun shifts the nstlist and rlist to
>> higher
>>> values. Isn't it likely to impact the simulation reliability in such
>> cases
>>> ?
>>>
>> Be careful with terms here. *Reliability* is whether the combination of
>> code and force field correctly reproduces the results for which it was
>> parameterized. Since you cannot use the twin-range implementation of the
>> group scheme in the Verlet scheme, you will be using a different setup. You
>> should absolutely evaluate reliability yourself - can you reproduce the
>> experimental observables that the parameterization claims, and do they
>> accompany plausible values for the observables in which you are interested?
>> The Verlet scheme uses and interprets rlist differently (and if we had our
>> decisions over again, we should not have been re-using rlist and nstlist
>> for schemes with rather different behaviour).
>>
>> The question of which setup is more likely to produce a useful physical
>> model for given computational cost is open, but consensus in the field
>> suggests that the features of
>>
>> a) symplectic integration (which implies reversibility),
>> b) suitable conservation of a quantity (thus adequate buffering),
>> c) explicit handling of non-isotropic long-range effects (including at
>> least electrostatic, but preferably also dispersion PME)
>>
>> are highly desirable (particulary if the system is not isotropic).
>> Unfortunately the scheme that Poger used in their parameterization has none
>> of those properties.
>>
>> Mark
>>
>>
>>> Thank you in advance for your replies!
>>> --
>>> 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.
>>>
>> --
>> 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.
>> --
>> 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.
>>

-- 
Dr Thomas Piggot
Visiting Fellow
University of Southampton, UK.