[gmx-users] ewald_rtol and convergence

[Mark Abraham]

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
> configuration.
 >
> Ewald_rtol                   total                  Short range       
>  Long range         1/beta (nm)
> 1e-02                  -53195.725708   -46946.948640     -6248.777068 
>      0.494129
> 1e-03                  -53432.220675   -45014.961538     -8417.259138 
>      0.386805
> 1e-04                  -53652.773007   -43093.541701   -10559.231306  
>     0.327146
> 1e-05                  -53872.777698   -41119.255332   -12753.522366  
>     0.288146
> 1e-06                  -54092.976378   -39070.390423   -15022.585955  
>     0.260198
> 
> 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
> energy?

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   
>    0.494129
> 1.8                  -52875.688467   -49910.307947     -2965.380520   
>    0.988258
> 2.7                  -52798.032604   -50836.767614     -1961.264991   
>    1.48239
> 3.6                  -52764.124490   -51296.887618     -1467.236872   
>    1.97652
> 4.5                  -52745.340510   -51572.825196     -1172.515314   
>    2.47065
> 
> Ewald: (fourier_nx,ny,nz = 25)
> 
> 0.9                  -52899.579997   -46946.948640    -5952.631357    
>     0.494129
> 1.8                  -52785.799065   -49910.307947    -2875.491118    
>     0.988258
> 2.7                  -52761.626727   -50836.767614    -1924.859113    
>     1.48239
> 3.6                  -52747.927053   -51296.887618    -1451.039435    
>     1.97652
> 4.5                  -52738.223387   -51572.825196    -1165.398191    
>     2.47065
> 
> 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?

Mark