[gmx-users] ewald_rtol and convergence
Mark.Abraham at anu.edu.au
Tue Feb 21 02:14:04 CET 2006
Michael Shirts wrote:
> Hi, all-
> So, I'm wondering what the effect of ewald_rtol is. My understanding
> is that it should just be shifting the calculation between the real
> space and the reciprocal space summations.
It also shifts the accuracy of the two components in a complementary
manner. See (for example) H.G. Petersen "Accuracy and efficiency of the
particle mesh Ewald method" J. Chem. Phys. 103 3668 (1995). Decreasing
fourierspacing and increasing pme_order (for PME) can improve the PME
reciprocal-space accuracy, as you know. There's an Essman J. Chem. Phys
paper worth looking at too.
> But the total
> electrostatic is changing non-negligibly as well. My system is a box
> of waters with a protein -- total charge is zero. PME Cutoff is 0.9,
> fourier_spacing is 0.02 (so well converged in all cases), PME order 6.
> Octahedral box, (5.78643,5.45550,4.72461) box vector, all the same
> Ewald_rtol total Short range
> Long range 1/beta (nm)
> 1e-02 -53195.725708 -46946.948640 -6248.777068
> 1e-03 -53432.220675 -45014.961538 -8417.259138
> 1e-04 -53652.773007 -43093.541701 -10559.231306
> 1e-05 -53872.777698 -41119.255332 -12753.522366
> 1e-06 -54092.976378 -39070.390423 -15022.585955
> I did some comparison to straight Ewald, as well.
> Ewald_rtol total Short range
> Long range
> 1e-02 -52899.579997 -46946.948640 -5952.631357
> 1e-03 -53001.090875 -45014.961538 -7986.129337
> 1e-05 -53182.993892 -41119.255332 -12063.738560
> We see that the short range part is exactly the same in PME and
> straight Ewald, as it should, though the long range part is not the
> same, even when the PME part is well converged (i.e., I reduced the
> fourier_spacing until any changes were sub 0.002 kJ/mol, which is one
> part in 25 million).
That's a convincing study of the phenomenon. It's worth remembering that
the accuracy of the individual components, their magnitudes, and the
values of rcoulomb, ewald_rtol, fourierspacing and pme_order are all
inter-related. In both of the above you are observing the accuracy
decrease (w.r.t. the presumed correct value of -52730) while varying
rtol. Since pme_order and fourierspacing are more than adequate, the
suspect parameter should be rcoulomb. I don't have an estimate for the
error in the real-space Ewald sum at hand, but the above reference cites
J. Kolafa and J.W. Perram, Mol Simul. 9 351 (1992) for an estimate of
the error in the reciprocal space *forces* proportional to
1/R * e^(-(beta * R)^2) for the cut-off R and Ewald parameter beta.
ewald_rtol * R = erfc (beta * R), so substituting a first-order
polynomial approximation for inverse erfc, we get the error proportional
to 1/R * e^(-(ewald_rtol * R)^2). This does increase with decreasing
ewald_rtol which is a bit weird though...
> I'm unsure why the differences. It there coding issue with the
> octahedral box? Is there some effect with the dielectic at infinity
> (i'm using epsilon_surface = 0). What is the 'true' electrostatic
I doubt there's a problem with the octahedral box, but if you wanted to,
you could test by using a cubic one. You can obviously approximate the
true electrostatic energy with a cut-off calculation with long cut-offs,
but that's roughly what you've done below anyway.
> I tried different cutoffs -- presumably, as we take a low rtol,
> pushing the calculation further into the direct space, and choose a
> longer cutoff, the correction will get more accurate.
> PME: (order 6, ewald_rtol = 1e-02)
> cutoff total Short range
> Long range 1/beta
> 0.9 -53195.725708 -46946.948640 -6248.777068
> 1.8 -52875.688467 -49910.307947 -2965.380520
> 2.7 -52798.032604 -50836.767614 -1961.264991
> 3.6 -52764.124490 -51296.887618 -1467.236872
> 4.5 -52745.340510 -51572.825196 -1172.515314
> Ewald: (fourier_nx,ny,nz = 25)
> 0.9 -52899.579997 -46946.948640 -5952.631357
> 1.8 -52785.799065 -49910.307947 -2875.491118
> 2.7 -52761.626727 -50836.767614 -1924.859113
> 3.6 -52747.927053 -51296.887618 -1451.039435
> 4.5 -52738.223387 -51572.825196 -1165.398191
> So it seems that -52730 is closer to the "right" answer. Why isn't
> this more possible to with other combinations of terms? Am I missing
> something? Is this level of agreement (both for Ewald, and between
> PME and Ewald) an unreasonable expectation? :)
That's a harder question to answer - it's a complex system. Maybe close
examination of those papers will shed some light?
More information about the gromacs.org_gmx-users