[gmx-developers] ewald_rtol / beta

[Mark Abraham]

On 26/06/2011 7:48 PM, Dommert Florian wrote:
> Hello,
>
> what about allowing to insert beta directly ?
>
> If a negative ewald_rtol is interpreted as the splitting parameter beta
> just two lines of codes have to be inserted after the declaration of the
> variables in calc_ewaldcoeff in ewald_util.c:
>
> if (dtol<0.0)
> return -dtol;
>
> What do you think about adding this feature to the code ?

Seems reasonable to me, but grompp should write a note to that effect, 
lest someone inadvertently use the feature.

Mark

> /Flo
>
>
> On Tue, 2011-06-07 at 11:28 +1000, Mark Abraham wrote:
>> 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
>> -- 
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-developers-request at gromacs.org.
>