Anistropic grids for Ewald Coulomb evaluation

[jtkrogel]

As noted in #2105, QMCPACK's periodic Coulomb energy calculation has huge errors for quasi-2D systems with large enough vacuum.

In #2125, a reference Ewald implementation was added to safeguard all runs against errors of this type via a reference check followed by an abort in the event of large errors.

The checking code was slow, requiring further optimization to be usable for large scale production runs. The code was brought to an acceptable speed across efforts contained in #2169, #2172, and #2176. The largest increases in speed were found by optimizing the gaussian exponents in the Ewald breakup, which minimized the number of points in the isotropic (cubic) sums.

While developing #2125, it became clear that many points within the cubic grids were contributing little for anisotropic quasi-2D systems. Efforts begun on adaptive anisotropic grids revealed that the number of anisotropic grid points required to sum the Ewald energy to high accuracy is on the order of the number of points already used by QMCPACK in production (which currently has a low accuracy), see #2182.

This raises the prospect that the main QMCPACK implementation of Ewald sums could be replaced with an adapted version of the anisotropic code. The resulting code would be much higher accuracy (resolving #2105) while likely remaining similar in speed (needs investigation).

The work plan here is to explore this possibility.

• 1. Investigate the size of anisotropic grids needed for full convergence of the Ion-Ion energy. Compare these to the underconverged spherical grids currently used by QMCPACK.
• 2. Explore the comparative speed of the anisotropic Ewald code to QMCPACK's optimized breakup in the context of the electron-electron interaction.
• 3. If higher accuracy can be obtained at comparable speed, investigate what effort would be required to fully replace QMCPACK's Ewald infrastructure with the anisotropic code.

[prckent]

Before we get too far in investigation or implementation, I think it would be prudent to discuss the goals of updating the Coulomb code and the conditions that should be met for the code to go in. This would apply to implementations and algorithms other than those in #2182 including any new ideas from the literature.

Goals

1. Compute all the Coulomb interactions (usually ion-ion, electron-ion, and electron-electron) to a specified tolerance, thereby being friendly to users and avoiding the recent surprises.

2. Be reasonably fast while meeting the accuracy criteria. e.g. No worse than 3x slower than the existing Optimized breakup code for dense materials where the Natoli-Ceperley implementation works well. Practically fast for the slow-convergence-with-old-algorithm cases identified so
   far such as graphene with 10 and 30 Angstrom gap.

3. Be completely unit testable with simple code. This will allow existing integration tests and initial Ewald check code to be largely done away with and replaced by faster and easier to modify and debug unit tests. They can also be more comprehensive. In practice this means that the guts of the implementation should minimize use of the existing QMCPACK class hierarchy.

4. Should be as simple as possible, avoiding premature optimization that reduces readability or that might not port to all architectures and compilers.

5. Replace the Optimized Breakup and non-reference Ewald codes. These should be deleted as the new code is proven. The new code will be more accurate and/or faster converging.

Since the initial testing of the adaptive anisotropic code suggests that performance is better for similar numbers of grid points and therfore numerical work, I am optimistic that all these goals can be met.

Others may wish to modify or add to the above goals.

I agree with Jaron's work plan that the next step for the anisotropic implementation is to do some timing tests to see where in the ballpark the current anisotropic code is performance wise compared to the production optimized breakup implementation.

[prckent]

Comments are still welcome, but since no one has suggested an alternative route, the next step "is to do some timing tests to see where in the ballpark the current anisotropic code is performance wise compared to the production optimized breakup implementation."