[gmx-users] confused about rcoulomb<=rlist

[Berk Hess]

>From: "David Mobley" <dmobley at gmail.com>
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: "Discussion list for GROMACS users" <gmx-users at gromacs.org>
>Subject: [gmx-users] confused about rcoulomb<=rlist
>Date: Wed, 12 Jul 2006 07:58:23 -0700
>
>Dear all,
>
>I desperately need some clarification about how GROMACS handles
>coulomb interactions with rcoulomb and rlist, and why it makes sense
>to do things this way.
>
>In particular, there was a thread on the revisions list
>(http://www.gromacs.org/pipermail/gmx-developers/2006-April/001581.html)
>a while back explaining that in 3.3.1, grompp checks to ensure that
>rcoulomb=rlist for PME/Ewald/PPM. I was very surprised by this, as I'd
>always been using rlist>rcoulomb, which I thought was best. Can
>someone explain why one would want rcoulomb=rlist? I don't quite
>understand the explanation on the revisions/developers list.
>
>My confusion is exacerbated by the fact that the GROMACS 3.3 manual
>includes the following statement:
>"The neighbor search cut-off rlist should be 0.1 to 0.3 nm larger than rvdw
>to accommodate for the size of charge groups and diffusion between neighbor
>list updates."
>
>This is exactly the same logic I was applying to rcoulomb when I
>concluded that I wanted to make sure rlist was larger than rcoulomb --
>basically, that I didn't want things outside rlist diffusing to within
>rcoulomb between neighbor list updates, and hence not being counted in
>the real-space interactions when they should be. Is there some other
>concern that outweighs this or means that this point is irrelevant?

I have added the check.

The problem was that Gromacs did not truncate the potential
at rcoulomb with PME, only rlist was used. However rcoulomb
was used in the determination of beta.
To avoid that people thought that rcoulomb had an effect
on the cut-off, I have added the check.

What is the optimal cut-off scheme is a different issue.
Indeed one would always want the force to go smoothly
at the cut-off, or before the cut-off in case one has charge groups
or nstlist>1.
However for PME one can not have 'exact' electrostatics while
the particle-particle force is zero at the cut-off since the reciprocal
space requires an error function contribution in real space.
Therefore the real space interaction has infinite range, but decays
very rapidly.
In PME the real space interaction is erfc(beta r)/r, which in Gromacs
is applied to all atom pairs in the neighborlist. The cut-off error
made here is very small, since ewald_rtol=1e-5 (I also often use 1e-6).
The alternative would be to somehow make the direct space interaction
go exactly to zero before the cut-off. This would lead to larger
errors, as one then misses more of the electrostatic interactions.
The only advantage would be that the integration could be more
reversible (assuming that all other algorithms are well reversible).

Given the accuracies of all other algorithms I would say it does not
make sense to remove some electrostatic interactions at the cut-off
when they are already 1e-5 or 1-e6 smaller than the Coulomb
interaction at that distance.

Berk.