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

[David van der Spoel]

Mark Abraham wrote:
>> 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.
It's center of geometry that is use as the criterion. For large charge 
groups this can be problem,, and therefore you want to add roughly the 
size of one charge group to the neighboursearching length.

> 
>> 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).
> 
I'd say avoid RF at all cost, see my recent paper and references 
therein: J. Chem. Theor. Comp. 2, 1-11 (2006), but for balance, also 
read some recent papers from the Hunenberger group.

-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++