Support for computing arbitrary order derivatives of 2-body integrals with Engine

[adabbott]

This PR adds support for 2-body integral Engines which compute 3rd and 4th derivatives. I have tested this current implementation up to d functions for standard 4-center operator::coulomb integrals only, and it is working as intended. Code has been added for 2-center and 3-center coulomb integrals as well, but not tested as my current libint configuration does not support their derivatives, and I do not have a reference implementation to test against. Hence, that part of the code is commented out for now.

@evaleev let me know if you want me to change anything. In particular,

• The lookup-arrays mapDerivIndexN_xxxx are very large and unsightly at higher orders, do you think they should instead be imported from another header file, and loaded in when the Engine object is instantiated with a particular deriv_order? Suggestions welcome
• I generated those new mapDerivIndexN_xxxx arrays with a short NumPy script, and this can be added to the PR if desired, or shared by other means.
• I am not sure how to implement the relevant unit tests

[evaleev]

@adabbott thank you, this looks like a great step towards the long-time goal!

Qs/comments:

• how hard would it be to port the logic of the script generating the maps to a C++ function? Ideally we could just do this for arbitrary orders ... the static maps were purely used to save time / indulge laziness.
• the map generation can be made purely static and triggered by libint2::initialize()
• ad hoc unit tests should go into tests/unit/test-2body.cc
• integration testing of libint is done via hartree-fock++, in the latter I would just compute the "core" contribution to the 2-body part of 3-rd/4-th deriv tensor and hardwire the reference results to hartree-fock++-validate.py

Thoughts?

[adabbott]

@evaleev Generating the maps with a C++ function sounds like a great idea. It shouldn't be too difficult to convert my python code to C++. I think I can write a universal function generate_deriv_index_map which given a deriv_order and BraKet type can output any particular mapDerivIndex* that currently appears in engine.impl.h. Do you have a suggestion as to where I would put this function? Also, where should the map array be created? Here, perhaps?

I'll look into adding tests once this is done.

[evaleev]

@adabbott great! I would make these tables managed by a dedicated singleton class whose entire purpose is to manage the initialization/lifetime and access from Engine ... see https://github.com/evaleev/libint/blob/master/include/libint2/solidharmonics.h#L82 for an example ... since Engine constructors can be called in general from multiple threads I would make the maps statically initialized up to the order that the library actually supports (this way you can avoid needing to lock the table of maps in case it needs to be extended to support derivatives of order higher than requested previously).

[adabbott]

@evaleev I think this is ready for another round of feedback. As this is my first substantial contribution to a C++ library, I am a bit out of my comfort zone, so nitpicking is very much welcomed. Some things that are worth noting:

1. Currently in initialize.h I have it calling DerivMapGenerator::initialize(), which computes and stores all of the mapDerivIndex arrays needed. This is working, but is this the correct place to put it?

2. Am I right in assuming that the definitions in engine.impl.h regarding LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri... can be removed and instead just always set int ket_lmax = hard_default_lmax_?

3. For DerivMapGenerator::instance I originally wanted to pass a BraKet type to this function, but including engine.h to give access to the enum class causes some strange errors. Currently I am instead passing an integer ncenters to distinguish between BraKet::xx_xx and BraKet::xs_xx, but this is probably not general enough as it does not distinguish between BraKet::xs_xx and BraKet::xx_xs, for example. I could change this function to accept a string, "xx_xx", "xs_xx", etc so that it is more easily extensible, but I'm open to better options. WDYT?

4. DerivMapGenerator::generate_deriv_index_map only generates the index mappings once when libint is initialized, for every BraKet type and deriv order up to LIBINT2_MAX_DERIV_ORDER. The overhead of doing this is minimal: generating all map arrays up to 6th-order derivatives for both Braket::xx_xx and BraKet::xs_xx takes about 300 ms and uses 4 kB on my machine.

[adabbott]

@evaleev I think this is ready for another review. All tests are passing locally when I configure --with-max-am=2 --with-opt-am=0 --enable-1body=2 --enable-eri=2 --enable-eri3=2 --enable-eri2=2, but for some reason the Travis build seems to be consistently failing when it gets to the pearl test scripts in src/bin/test_eri/ ... any idea why?

[evaleev]

the failing CI job in https://github.com/evaleev/libint/actions/runs/1233407450 fails due to GH actions jitter, rest looks good.