[gmx-users] The meaning of rlist and rcoulomb for coulombtype=reaction-field

[Mark Abraham]

> I am having a difficult time understanding the meaning of 'rlist' and
> 'rcoulomb' when the reaction-field treatment of electrostatics is used.
>
> After explaining the twin-range cutoff method (which seems to be
> fairly straightforward) which explains the case where rlist <
> rcoulomb, the manual states:
>
> "Except for the plain cutoff, all of the interaction functions in
> Table 4.2 require that neighbor searching is done with a larger
> radius than the r_c specified for the functional form, because of the
> use of charge groups. The extra radius is typically of the order of
> 0.25 nm"
>
> This appears to suggest that for reaction-field, we would wish to set
> rlist = rcoulomb + 0.25 nm, something which grompp does not allow.

7.3.8 specifies rlist as the short-range neighbour list, so you definitely
don't want your suggested equation to hold. As written, I think that
exerpt implies that you want to use rlist+delta and rcoulomb+delta for
some delta for the neighbour searching. However that doesn't make complete
sense in the context of charge groups, because you have to decide what to
do with charge groups that straddle the new boundary as well...

>  From my experience with molecular simulation, it seems to be
> desirable to ensure that the neighbor list maintained between the
> times neighbor searching is conducted is indeed larger than the
> cutoff beyond which electrostatic interactions are zero.

That's correct for single-range, where the neighbour list is used to prune
the search space to find the actual atoms inside the cut-off for an
energy/force evaluation, I believe.

With twin-range, as described in the manual just before your exerpt quoted
above, periodically you have the short-range neighbour list constructed
from those pairs inside rlist (plus perhaps charge groups extending
outside), and if max(rcoulomb,rvdw)>rlist then there's an energy and force
evaluation for that near-spherical shell between rlist and
max(rcoulomb,rvdw) that is evaluated then and applied to each integration
step before the next neighbour list construction. Between constructions,
the force and energy contributions of the inner sphere are evaluated each
time. The working assumption here is that neighbour list construction is
frequent enough that you don't have substantial distortion of that inner
sphere between constructions. This is a quite different use of the
neighbour list, IMO, since there is no "pruning" stage.

> Does this 0.25 nm buffer get automagically added to the specified
> rcoulomb to determine the neighbor searching radius?

No, I checked, and that doesn't happen in src/mdlib/ns.c

My guess is that the actual implementation includes all atoms in charge
groups that have at least one atom inside the respective rlist and
rcoulomb cut-offs, and that accordingly the documentation is out of line,
but I don't care enough to read ns.c that closely to find out, so
hopefully David/Berk/Erik can clarify here.

> Also, what
> happens when rlist < rcoulomb for reaction-field?  Is the standard
> twin-range method employed here?

Yep, that's how it tells whether twin-range or simple is desired.

> Finally, I have very little
> experience with reaction-field -- is there a recommended combination
> for rlist and rcoulomb for simulations of small protein systems?

Heh... thorny issue. One way to maximise the applicability of your force
field to your problem is to use the same cut-off regime that it did when
it was being parameterized. However, some older force fields used cut-offs
imposed by available computing time that are significantly smaller than
you could use on today's computers, so it's probably defensible to extend
the cut-offs. This will have a small and unknown perturbation on how well
your force field now models real physics... it hopefully has the right
sign... but if not, that problem is probably swamped by the general
unsavouriness of the point-charge assumption in the first place, anyway.

I'd suggest looking around at recent papers that use the force field
you're working with, on the sort of size problems you're working on, with
particular focus on the papers by the principal authors of that force
field (since they can know things about the force field that haven't made
the literature).

Mark