[gmx-developers] ewald_rtol / beta
Mark Abraham
Mark.Abraham at anu.edu.au
Tue Jun 7 03:28:30 CEST 2011
On 7/06/2011 7:55 AM, Dimitris Dellis wrote:
> Hi gmx_developers.
>
> In manual (v4.5) p. 100 ewald_rtol parameter is referred as "The
> ewald_rtol parameter is the relative strength of the electrostatic
> interaction at the cut-off"
> According to this (and some older posts in lists), beta parameter is
> obtained by solving equation
> ewald_rtol = erfc( beta * r_c) / r_c [1]
> In gmxlib/ewald_utils.c, function calc_ewaldcoeff() (versions 3.3.4 -
> 4.5.4) beta is obtained by solving eq.
> ewald_rtol = erfc( beta * r_c) [2]
> (division by r_c omitted ).
> These 2 expressions give the same beta only if r_c=1.
Yep - though the division by r_c in [1] is erroneous. Apologies - I was
probably one of the posters to whom you refer :-)
If ewald_rol is the "relative strength of the electrostatic interaction
at the cut-off" then we have
ewald_rtol = V_Ewald(r_c) / V_normal(r_c)
so
ewald_rtol = [erfc(beta * r_c ) / r_c ] / [1 / r_c]
= erfc(beta * r_c )
In order to determine beta, one could equally well define the .mdp
parameter ewald_maxv to satisfy
ewald_maxv = erfc(beta * r_c) / r_c
or even have beta as an .mdp parameter. All that is necessary is that
some algorithm exists to uniquely define beta, perhaps as a function of
r_c - it doesn't matter what it is.
> I tried to investigate convergence of real/recip. sums as function of
> beta. For certain (beta, r_c), I calculate ewald_rtol using [1], put
> this ewald_rtol to mdp.
> In md.log the reported beta is not the same with the one that was used
> to calculate ewald_rtol (not only last few digits deviation), except
> if r_c=1.
> If I calculate ewald_rtol using [2] and give this value in mdp, mdrun
> reports the expected beta.
Yep.
> In almost all systems this doesn't affect results since real/recip.
> sums converge for values of beta slightly different (few %) than the
> supposed beta,
> and typical rcoulomb is close to 1.
> This might affect results in some cases that exhibit difficult (if not
> impossible) real/recip. sums convergence, that for various reasons
> require large cut-offs.
>
> Is this a bug ?
I think only in the sense that incomplete documentation leads to
incomplete understanding. Even if the original coder intended ewald_rtol
to have the sense of ewald_maxv, the resulting code serves the purpose
in a slightly different way. You can see in src/mdlib/tables.c that
ewaldcoeff is used in the erfc(ewaldcoeff * r)/r sense, as expected.
Assuming nobody corrects me, I'll stick a comment in the code and a
clarifying sentence in the manual.
Mark
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20110607/187ae23f/attachment.html>
More information about the gromacs.org_gmx-developers
mailing list