[gmx-users] Gromos54a7 electrostatics interactions ?

[Piggot T.]

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.