[gmx-developers] Re: [gmx-users] confused about rcoulomb<=rlist

[David Mobley]

Berk,

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

OK, so, in other words, rcoulomb now only affects beta with PME; rlist
is used for the electrostatic cutoff.

Is there a straightforward way this could be changed in future
versions so that the cutoff used for electrostatics could be different
from rlist? In particular, I am concerned that this means that one
can't vary the vdw cutoff independently of the electrostatics cutoff
beyond some limit. That is, rlist needs to be >rvdw with shift/switch
vdw potentials:

WARNING 1 [file equil_constp0.7.mdp, line unknown]:
  For energy conservation with switch/shift potentials, rlist should be 0.1
  to 0.3 nm larger than rcoulomb/rvdw.

Suppose I want to run with a 1.4 nm rvdw -- then I now have to use
rlist=1.5 nm and rcoulomb=1.5 nm? That doesn't seem good. One would
like to be able to change rvdw independently of the electrostatics --
for example, to explore the effects of long range van der Waals
cutoffs. In fact, Michael Shirts and I and others had been trying to
do exactly that in GROMACS 3.1.4 and 3.3 (adjust rvdw and rlist
independently of rcoulomb in binding free energy calculations, in
order to see how strongly the choice of vdW cutoff affects computed
binding free energies). If changing rlist also changes the real space
cutoff for electrostatics, things are substantially more complicated
and it seems unclear how to even address this issue.

Maybe a solution would be to NOT use rlist for beta and for the PME
potential, and instead to use rcoulomb or some extra variable... Then
one could change the neighbor list (rlist) and rvdw without having to
tweak the electrostatic interactions...

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

So, if I understand properly, you're saying that the notion of an
electrostatics "cutoff" doesn't really make sense with PME, so my
logic about having the neighborlist extend past the "cutoff" doesn't
really apply. Rather, one just wants to ensure that the "cutoff" is
sufficiently large that electrostatic interactions are fairly small
and thus any inaccuracies due to neighbor list issues will be
basically negligible...?

Thanks,
David


> Berk.
>
>
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>