[gmx-developers] ewald_rtol / beta

[Mark Abraham]

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>