EKinetic Gradient from LibInt

[ryanmrichard]

This issue depends on NWChemEx/SimDE#40 and https://github.com/NWChemEx-Project/Mokup/issues/31.

Assuming the description in NWChemEx/SimDE#40 is correct, Libint will give us quantities like (one shell-pair at a time). This issue focuses on writing a module LibintEKineticNuc which will:

• Satisfy the EKinetic_Nuc PT defined in Property Type for Gradient of EKinetic SimDE#40
• Initialize Libint (see the tutorial link at the bottom of this issue)
• Create a lambda function which can be used to fill in a TA::DistArrayD instance via TA::make_array (the signature for the lambda is double (TA::Tensor<double>& t, const TA::Range& range) where the return is the norm of the tile, t is the tile being filled in, and range is the range of the tile (the indices of the elements in t).
• Copy the TA tensor into a TW tensor and return the TW tensor.

Depending on when this issue is being worked on, NWChemEx/TensorWrapper#33 may have been resolved. If that's the case the lambda signature will change slightly and the copying of the TA tensor into the TW tensor will no longer need to happen.

The literal tasks:

• Write the LibintEKineticNuc module (I recommend putting it in a new file src/integrals/libint/ekinetic_nuc.hpp)
• unit test the module (I recommend putting the test in a new file tests/integrals/libint/test_ekinetic_nuc.cpp) the testing should be largely analogous to how we unit test EKinetic itself (see here)

Additional Notes:

• We have some infrastructure in this repo for initializing and calling Libint. That infrastructure is heavily templated, and frankly, at this point a mess (it's built off a legacy version of what's now called PluginPlay and can be massively streamlined). It's recommended you ignore most of that infrastructure and follow the tutorials for Libint to implement LibintEKineticNuc
• For the purposes of this issue there's no need to worry about attempting to generalize the module for different operators, just hard-code everything to the electronic kinetic energy. We'll worry about generalizing in another issue.
• While I recommended ignoring most of the infrastructure in this repo, some of it is likely to still be very useful, particularly take a look at nwx_libint.hpp which will help convert from Chemist objects to Libint.
• As a first pass just worry about an implementation where the entire integral is stored in core memory. Direct is done by changing the tile types of the TA tensor and it should be possible to implement direct in a manner that is largely orthogonal to the filling in of the tensor handled by this module.
• At this stage I wouldn't worry too much about documentation. Any documentation that is written would probably need changed when the module gets generalized to other derivative types.

[hjjvandam]

According to https://github.com/evaleev/libint/wiki/using-modern-CPlusPlus-API the key to setting up the right integrals is in creating the right integral Engine. The Engine constructor can take an optional argument for the order of differentiation with respect to the nuclear coordinates.
For integral derivatives the constructor of the integral engine has the following API:

Engine(Operator oper, size_t max_nprim, int max_l, int deriv_order = 0);

Here

• oper is the value of an enumerator of operators,
• max_nprim is the maximum number of primitives in the contraction of Gaussians in the basis set,
• max_l is the maximum angular momentum in the basis set,
• deriv_order is the order of differentiation, by default it is 0, but it can be set to an arbitrarily high integer number.

The operators available and relevant to this issue are:

• Operator::overlap for overlap integrals
• Operator::kinetic for kinetic energy integrals
• Operator::nuclear for nuclear attraction integrals
• Operator::coulomb for electron repulsion integrals

[hjjvandam]

Note that by default Libint supports derivatives of at most order 0. Calculating integral derivatives of higher orders (like 1) requires regenerating Libint.

[hjjvandam]

Libint generates integral derivatives for a given set of shells. To store these derivatives in the right place we need to know the ranks of the corresponding atoms in the overall geometry. In Libint the function

    std::vector<long> shell2atom(const std::vector<Atom>& atoms)
    static std::vector<long> shell2atom(const std::vector<Shell>& shells, const std::vector<Atom>& atoms, bool throw_if_no_match = false)

The result shell2atom[ishell] the rank of the atom that shell ishell resides on.

[hjjvandam]

To access the integrals we need access to the memory where they are stored:

     const auto& Engine.results()

This points to a linear array in row-major form. I.e. int(f1,f2) is stored in results()[f1*n2+f2].
The evaluation of the integrals is performed by

    Engine.compute(shell1,shell2,...)

Also note that to ensure all Cartesian basis functions are normalized to 1.0, one needs to call

     Engine::set(CartesianShellNormalization::uniform)

Derivatives are stored as 6 sets corresponding to d/dx1, d/dy1, d/dz1, d/dx2, d/dy2, d/dz2 for each derivative all the corresponding integral derivatives are stored. I.e. Engine.results() points to an array of pointers where Engine.results()[3] points to an array with the d/dx2 derivative integrals (in the case where the operator does not depend on the atomic coordinates).