[gmx-developers] ewald_rtol / beta

[Dommert Florian]

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 ?

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


-- 
Florian Dommert
Dipl. - Phys.

Institute for Computational Physics
University Stuttgart

Pfaffenwaldring 27
70569 Stuttgart

EMail: dommert at icp.uni-stuttgart.de
Homepage: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert

Tel.: +49 - (0)711 - 68563613
Fax.: +49 - (0)711 - 68563658
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 198 bytes
Desc: This is a digitally signed message part
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20110626/de5efd47/attachment.sig>