[gmx-users] Gromos54a7 electrostatics interactions ?

Mark Abraham mark.j.abraham at gmail.com
Fri Sep 22 13:41:49 CEST 2017


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

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.


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

More information about the gromacs.org_gmx-users mailing list