[WIP] Add cusp correction based on atomic orbitals
This implements the cusp correction scheme of Manten and Lüchow ( J. Chem. Phys. 115, 5362–5366 (2001) https://doi.org/10.1063/1.1394757 ).
The existing cusp correction scheme in QMCPACK is based on correcting molecular orbitals. This does not interact well with orbital rotation (one would need to recompute the correction after every rotation). Correcting individual atomic basis functions would be a better fit with orbital rotation.
This PR uses a script ("compute_cusp_correction.py") to compute the parameters and write them to an HDF file. The QMCPACK code reads those parameters and evaluates corrected orbitals.
This is a minimal implementation to prototype the idea.
The design should be improved for actual use. The GaussianCombo class should not be modified. Most likely a wrapper class (GaussianComboWithCusp?) should be created that evaluates the cusp function for r < r_c, and calls GaussianCombo otherwise. The one tricky part is that not all the orbitals needs to modified - only the 1s- and 2s-type orbitals. The storage of cusp parameters in the HDF file may need to change to reflect the design.
It has been tested with C_2 and cc-pv{d,t,q,5}z basis sets, and works with orbital rotation.
Testing with H2O yielded some problems:
- H - the second derivative of the cusp doesn't intersect the second derivative the GTO unless r_c is enlarged to 3.0-4.0. The paper uses molecules containing hydrogen, but makes no mention of special handling for it.
- O - H2O showed a larger variance than expected, so I tried an O atom. The O atom causes trouble for the optimizer - the energy and variance decrease for a few steps, then start increasing (sometimes dramatically). The problem is more severe with J3 included, but can be reproduced with just J1 and J2. This problem can also occur with the MO cusp, so it is not exclusive to the cusp correction scheme in this PR. Increasing minwalkers to larger value (0.8) seems to help, but this needs more investigation. (Note there has been no orbital rotation involved here - just optimizing the Jastrows)
What type(s) of changes does this code introduce?
Delete the items that do not apply
- New feature
Does this introduce a breaking change?
- No
What systems has this change been tested on?
desktop
Checklist
Update the following with a yes where the items apply. If you're unsure about any of them, don't hesitate to ask. This is simply a reminder of what we are going to look for before merging your code.
- Yes/No. This PR is up to date with current the current state of 'develop'
- Yes/No. Code added or changed in the PR has been clang-formatted
- Yes/No. This PR adds tests to cover any new code, or to catch a bug that is being fixed
- Yes/No. Documentation has been added (if appropriate)