adding qmcfinitesize tool

[camelto2]

Proposed changes

This PR adds the qmcfinitesize tool that implements some of the finite-size corrections from Holzmann et al.. This has been described in the manual but has never been in the main code. I updated the code from @rcclay to be up to date to develop. This requires some postprocessing of the S(k) from the SkAll estimator in order to use (for example, twist averaging and creating the extrapolated estimator instead of just using the mixed estimator).

This would close #1007 and #2000.

A following PR will be given which will implement Eqn. (35) to give the finite-size correction for the kinetic energy, provided a kSpace jastrow is included to capture the long-range behavior.

What type(s) of changes does this code introduce?

• New feature
• Code style update (formatting, renaming)

Does this introduce a breaking change?

• No

What systems has this change been tested on?

my workstation

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. This PR is up to date with current the current state of 'develop'
• Yes. Code added or changed in the PR has been clang-formatted
• No. This PR adds tests to cover any new code, or to catch a bug that is being fixed
• No. Documentation has been added (if appropriate)

[rcclay]

OK to test.

[qmc-robot]

Can one of the admins verify this patch?

[Paul-St-Young]

@camelto2 @rcclay Great to see this being merged in!

Have you guys tried using this version in a non-cubic box?
Back when I tried an old version, the 3D spline of S(k) would evaluate to garbage at large k (without throwing an error).
I had to guard eval_UBspline_3d_d by hand. See this commit.

I also found it beneficial to allow different long-range breakup schemes (Ewald, Natoli, and Esler), as well as different methods to calculate the thermodynamic value:

1. Integral using spherically average S(k) (current implementation)
2. 1/N extrapolation over 2 box sizes by directly using the 3D spline of S(k)
   They can serve as a robustness check/error estimate of the finite-size correction.

What are your future plans with this tool?

[camelto2]

@Paul-St-Young I have checked a few non-cubic cells and didn't notice any pathologies. However, I will look into adding some of the guards you mention in case they do come up.

As for the breakup, the choice of LR breakup schemes here is determined from the input xml file. So by default it is esler. You can change that by adding into your particle set
ewald/opt_breakup/opt_breakup_original
for ewald, esler, natoli respectively.

Currently we only have the one spherical average integral. I don't think there are plans to incorporate others for the time being. However, I will be adding some functionality to get errorbars on the correction by resampling the S(k) from the error bars given. If we determine other implementations are necessary we will add them later

[rcclay]

@Paul-St-Young For selecting the long range handler, the finite-size tool inherits functionality enabled by PR #2190 .

[rcclay]

@Paul-St-Young Excluding numerical issues with fitting the 3d splines... on a boundary for example... did you ever find a case where a radially averaged then integrated S(k) gave "better" results than the 3D spline integration?

[Paul-St-Young]

@camelto2 great the non-cubic cells are OK. Maybe the 3D spline has been updated since I tried this. In any case, the extra guards shouldn't do any harm.
It would be great if bootstrap error estimate became an option.

Regarding the choice of long-range handler, #2190 addressed it nicely.

@rcclay I never encountered a case where one approach was conclusively "better" than the other. Mostly, the two approaches agree with each other within say 0.1 mha/atom or less when neither runs into numerical issues. I personally prefer doing two boxes then extrapolate using 1/N, because this approach does not use splines so much and does not get too close to k=0, where numerical instabilities lurk.
Although, with bootstrap resampling error estimates, one should be able to spot instability of the splines, so this will (hopefully) become a non-issue.

[Paul-St-Young]

@rcclay quick question: after #2190, how can I use Esler break for total energy and Natoli break for the forces?

[rcclay]

@Paul-St-Young You don't anymore, sadly. The forces/energies use the same handler.

[prckent]

@Paul-St-Young We are quite minded to remove the "optimized" break up methods if we can implement a sufficiently fast and robust version of conventional Ewald. (Not guaranteed, but WIP #2182 may be that implementation). Several people were badly burned by the surprising poor convergence behavior of the "optimized" approaches for certain geometries, which basically makes them unusable as a default in a general code. Is there a strong reason to keep the other approaches around?

[rcclay]

Going to prevent merge until Cody works through some remaining issues/enhancements.

[ye-luo]

Is it easy to create some tests to monitor the health of added codes?

[prckent]

@ye-luo I second this question.

[camelto2]

@rcclay I refactored your old code quite a bit and tried to remove all of the unused code, but I may have missed a few things. I rearranged the correction calculations so they are packaged into smaller chunks, and that allowed me to resample the splines based on the errorbars provided in the input S(k) file. This code now creates error estimates for the various corrections based on these resampled splines. You may want to have a look at the changes, as things have been modified quite a bit.

@ye-luo @prckent Some sort of testing can be cooked up. It may take some time to get to it though

[ye-luo]

src/QMCTools/QMCFiniteSize/CMakeLists.txt seems like an empty file. Please remove it if so.

[ye-luo]

Could you check if the dependencies here are necessary?
Do you need qmcwfs or just einspline?
Is qmc target really necessary?

[camelto2]

at the moment, it only needs einspline. this will change in the future when a kinetic energy FS correction is added, which will depend on the 2body jastrow and a kspace jastrow. if it is preferable to just add the various .cpp files in FSSRCS that each thing depends on I can do that instead as opposed to targetting libraries that include those dependencies

[ye-luo]

Let us put einspline and change it to qmcwfs when actual WFS dependence comes in.

[rcclay]

For user friendliness, this tool should work with any QMCPACK xml input file. Unfortunately, even if you're not using any wavefunction specific information, the system initialization path that allows you to grab a particle set and box from ESHDF file, in the event it's not specified in the XML, needs to go through the wavefunction builder. So I think that qmcwfs is a necessary dependency right now.

[ye-luo]

Thanks. I see. The dependency is needed right now.

[ye-luo]

public QMCAppBase I'm not convinced that this inheritance is necessary. Do you use any thing from QMCAppBase class? Or it is just for convenience?

[camelto2]

we are parsing the input file to the qmcpack run that generated the S(k) to do the correction. We do that right to get particleset info and eventually the wavefunction stuff. Those things are defined inside of QMCAppBase, which is why it was derived from it

[ye-luo]

Does it parse full qmcpack xml input? If yes, it makes sense.

[camelto2]

it parses the xml to determine the simulation cell, generate a particleset, and wave functions. it doesn't care about the method sections, but basically everything else is parsed. After discussing this with @rcclay, qmcwfs is likely a hard requirement. In some xml inputs, one may just specify an hdf5 file without any of the cell information as part of the xml file. In that case, stuff is allocated by the wfnbuilder.

[camelto2]

updates to the tool. I have added unit tests, and compare against an analytic S(k) to make sure we get the right answer.

The unit tests add a directory in QMCTools so that anything in QMCTools currently or added in the future can easily be unit tested.

An extra function had to be added to LRHandlers to evaluate v_lr(k) at arbitrary k values, so that we would have an improved integrand for the integral part of the correction. Thanks to discussions with @rcclay and @Paul-St-Young for pointing it out.

[rcclay]

Test this please.