Reference Ewald implementation

[jtkrogel]

This PR adds a reference Ewald implementation with the intent to check the accuracy of QMCPACK's periodic ion-ion energy at runtime.

The reference Ewald implementation is complete and correct for the cases I've tested. I've checked it against QE and related Python implementations for VO2, graphite, and graphene (large c-axis cell).

The PR is marked as WIP to discuss what else should be done before adding a hard abort if QMCPACK and the reference disagree and merging into mainline. Possibly more extensive tests of the reference code should be performed to avoid false positives (i.e. it is bad form to block a user's run due to faulty reference behavior).

[qmc-robot]

Can one of the admins verify this patch?

[jtkrogel]

QMCPACK vs. reference for rutile VO2, graphite and graphene:

rutile VO2

  Reference ion-ion energy: -148.398013749
  QMCPACK   ion-ion energy: -148.39804090736
            ion-ion diff  : -2.7158357454482e-05

graphite

  Reference ion-ion energy: -15.770484019514
  QMCPACK   ion-ion energy: -15.77048578662
            ion-ion diff  : -1.767106002859e-06

graphene

  Reference ion-ion energy: 51.653202149774
  QMCPACK   ion-ion energy: 51.451335009833
            ion-ion diff  : -0.20186713994084

[rcclay]

OK to test

[prckent]

Thanks. I will have some comments in the week after Thanksgiving. Will help with #2105 .

[jtkrogel]

After discussion with @prckent, I added a hard abort when a difference of about 0.01 eV/atom is exceeded between QMCPACK's and the reference ion-ion energy. When the tolerance is exceeded, the user receives these instructions:

Error in ion-ion Ewald energy exceeds 0.0003 Ha/atom tolerance.

  Reference ion-ion energy: 51.653202149942
  QMCPACK   ion-ion energy: 51.451335009471
            ion-ion diff  : -0.20186714047129
            diff/atom     : -0.10093357023564
            tolerance     : 0.0003

Please try increasing the LR_dim_cutoff parameter in the <simulationcell/>
input.  Alternatively, the tolerance can be increased by setting the
LR_tol parameter in <simulationcell/> to a value greater than 0.0003. 
If you increase the tolerance, please perform careful checks of energy
differences to ensure this error is controlled for your application.

Fatal Error. Aborting at ion-ion check failed

I exposed the tolerance as an input parameter in case the user has a good reason to increase it. Later on, the tolerance can be used to fully replace LR_dim_cutoff as the means to control the LR sums.

[jtkrogel]

Ready for review.

[jtkrogel]

The CI is failing at a deterministic test of the PBCAA energy for a simple cubic Hydrogen cell (see test_coulomb_pbcAA.cpp).

The test for the H-H interaction passes:

REQUIRE(val == Approx(-1.418648723)); // not validated

While the test for the e-e interaction fails:

REQUIRE(val == Approx(-1.3680247006)); // not validated

Error in ion-ion Ewald energy exceeds 0.0003 Ha/atom tolerance.

  Reference ion-ion energy: -1.4186487397403
  QMCPACK   ion-ion energy: -1.368024558717
            ion-ion diff  : 0.050624181023291
            diff/atom     : 0.050624181023291
            tolerance     : 0.0003

In the H-H and e-e tests, the cell is the same and of course the species differ only by the sign of the charge, and so I would expect the Ewald energy to be the same (-1.4186487397403 Ha).

The test for the e-e part has this comment:

// Self-energy correction, no background charge for e-e interaction

Does anyone know if a term is being intentionally excluded here?

[jtkrogel]

@rcclay @markdewing Any thoughts?

[rcclay]

OK. Knee jerk reaction is to agree with you that the terms should be the same, but I'm going to have to carefully work through what the code is actually doing to be conclusive on that. Random comment on PR: please document src/QMCHamiltonians/EwaldRef.h, as that's the substantive piece of code in this PR. Like, what are the introduced auxiliary classes doing, what are the internal variables, what are the function calls doing specifically, etc.

[jtkrogel]

@rcclay, sufficient?

[rcclay]

Much better. Thank you :)

[markdewing]

About the comment "no background charge for e-e systems". In the python script (in qmc_algorithms/LongRange/ewald_sum.py) a comment mentions that the background charge is only non-zero for charged systems. An isolated e-e system would need the background charge to be correct. However, I suspect this term shouldn't be needed when all the terms are added together (e-e, e-n, n-n)??

Without the background charge, the values depend on the splitting parameter. The python script implements the Ewald algorithm, so values in test_coulomb_pbcAA_ewald.cpp can be reproduced with it.
The following changes to the script

1. The splitting parameter (kappa) is set to 5.47723 (found in the output when running the unit test)
2. Electrons only (qs=[-1.0], rs = [re1])
3. Set bg=0
   Then the script gets close to the values in the unit test (self energy correction = 3.090196 and the sum is -1.3662889)

We don't have a independent implementation of the optimized break-up method, so the values in test_coulomb_pbcAA.cpp will be difficult to reproduce.

[rcclay]

@markdewing Wait, is that what the code does--keep the BG term for ion-ion but omit it for e-e?

[markdewing]

@rcclay I don't know about the code. The analysis came from comparing the python script with the unit tests. It's still mysterious as to why the n-n result seems to include the background term, but e-e result does not.

[rcclay]

@markdewing That's concerning. I'll do some code spelunking.

[markdewing]

The difference in H-H and e-e comes from CoulombPBCAA::evalConsts, but I don't know where. My guess is that it is something to do with the set up of the SpeciesSet in the unit test.

[prckent]

Private class member data please.

[prckent]

Thanks for the update. Of course the unit test and ion-ion vs e-e issue will need solving before merge.

[jtkrogel]

@markdewing The values you find with the Python script agree much more closely with the test values used for the Ewald implementation (see test_coulomb_pbcAA_ewald.cpp):

  // Self-energy correction, no background charge for e-e interaction
  double consts = caa.evalConsts();
  REQUIRE(consts == Approx(-3.090194));

  double val = caa.evaluate(elec);
  REQUIRE(val == Approx(-1.366567)); // not validated

I expect the exclusion of the background term is what is causing the unit (regression) test to fail.

It actually seems like the issue uncovered here is separate from this PR. However, it probably motivates using the reference code to check more than just the ion-ion part during initialization. Maybe we should also check the total (e-e + e-I + I-I) part at the first Hamiltonian evaluation?

[markdewing]

The issue with the e-e unit test is the lack of setting the "membersize" attribute.

Adding to the following to TEST_CASE("Coulomb PBC A-A elec Ewald3D", "[hamiltonian]") in test_coulomb_pbcAA_ewald.cpp gives the same answers as the ion-ion test.

  int pMembersizeIdx           = tspecies.addAttribute("membersize");
  tspecies(pMembersizeIdx, upIdx) = 1;
  tspecies(pMembersizeIdx, downIdx) = 0;

[jtkrogel]

OK. Do you want me to update the unit tests this way?

[markdewing]

Yes. I think the same change should fix the non-Ewald tests as well.

[markdewing]

To be clear about the unit test issue - if the species membersize is not set, it ends up as zero, which sets NofSpecies to zero in CoulombPBCAA.cpp::evalConsts. Consequently, this makes the contribution of the background term zero (last loop in evalConsts). With the test fixed, the values in the unit test should be independent of the splitting parameter (to zero order. There's still convergence, which the reference implementation is checking).

[jtkrogel]

How do you get ctest to print the header of the test it is currently executing?

TEST_CASE("Coulomb PBC A-A elec", "[hamiltonian]")

I tried with --verbose but did not see this information.

[markdewing]

The "-s" option, which also prints successful tests will print the test name, but it might follow the test output.

[markdewing]

In case anyone doesn't know this, unit tests can be run individually by putting the test name on the command line. This is extremely useful when focusing on a particular test.

Most the names have spaces, so the name needs to be put in quotes.
tests/bin/test_hamiltonian "Coulomb PBC A-A elec Ewald3D"

[ye-luo]

In case anyone doesn't know this, unit tests can be run individually by putting the test name on the command line. This is extremely useful when focusing on a particular test.

Most the names have spaces, so the name needs to be put in quotes.
tests/bin/test_hamiltonian "Coulomb PBC A-A elec Ewald3D"

@markdewing could you document this in our manual? I didn't know this feature and always went to modify CMakeList to reduce the number of tests for debugging. Thank you.

[markdewing]

It's mentioned in the manual in the chapter on unit tests. Though perhaps it could use an example to make it stand out more.

[jtkrogel]

There were a number of bugs of that type in the unit tests. All should be fixed now.

[prckent]

@markdewing Example in the manual would help. It will make the result more discoverable through search.

[prckent]

I will merge this now that the weekly integration tests have run.