Libint2 gives some two-electron integrals that are permutation-symmetry-inconsistent?

Hello!

Thank you for this amazing electron integral library!

I've recently encountered a strange behaviour when trying to compute the two-electron integrals for a stripped-down STO-3G basis of a single vanadium atom that consists of only one P shell (the first P shell in STO-3G) and one D shell (the only D shell in STO-3G), both of which are set to uniformly normalised Cartesian Gaussian functions:

V     0
P   3   1.00
      0.8378385011D+02       0.1559162750D+00
      0.1946956493D+02       0.6076837186D+00
      0.6332106784D+01       0.3919573931D+00
D   3   1.00
      0.2964817927D+01       0.2197679508D+00
      0.9043639676D+00       0.6555473627D+00
      0.3489317337D+00       0.2865732590D+00
****

The code for this is

#include <libint2.hpp>
#include <stdlib.h>

using libint2::Shell;
using libint2::BasisSet;
using libint2::Operator;

int main(int argc, char* argv[]) {

    std::vector<libint2::Shell> shells;
    shells.push_back(
        {
            {0.8378385011e+02, 0.1946956493e+02, 0.6332106784e+01}, // exponents
            {   // P shell, spherical=false, contraction coefficients
                {1, false, {0.1559162750, 0.6076837186, 0.3919573931}}
            },
            {{0, 0, 0}} // origin coordinates
        }
    );
    shells.push_back(
        {
            {0.2964817927e+01, 0.9043639676e+00, 0.3489317337e+00}, // exponents
            {   // D shell, spherical=false, contraction coefficients
                {2, false, {0.2197679508, 0.6555473627, 0.2865732590}}
            },
            {{0, 0, 0}} // origin coordinates
        }
    );

    libint2::BasisSet obs = libint2::BasisSet();
    for (libint2::Shell shell : shells) {
        obs.push_back(std::move(shell));
    }
    auto shell2bf = BasisSet::compute_shell2bf(obs);

    libint2::initialize();
    libint2::Engine engine = libint2::Engine(libint2::Operator::coulomb,
                                             libint2::max_nprim(obs),
                                             libint2::max_l(obs),
                                             0);
    // Force uniform normalised Cartesian functions
    engine.set(libint2::CartesianShellNormalization::uniform);
    const auto& buf = engine.results();

    for(auto s1=0; s1 != obs.size(); ++s1) {
        const auto bf1_first = shell2bf[s1];  // first basis function in this shell
        auto n1 = obs[s1].size(); // number of basis functions in shell 1
        for(auto s2=0; s2 != obs.size(); ++s2) {
            const auto bf2_first = shell2bf[s2];  // first basis function in this shell
            auto n2 = obs[s2].size(); // number of basis functions in shell 2
            for(auto s3=0; s3 != obs.size(); ++s3) {
                const auto bf3_first = shell2bf[s3];  // first basis function in this shell
                auto n3 = obs[s3].size(); // number of basis functions in shell 3
                for(auto s4=0; s4 != obs.size(); ++s4) {
                    const auto bf4_first = shell2bf[s4];  // first basis function in this shell
                    auto n4 = obs[s4].size(); // number of basis functions in shell 4
                    engine.compute2<Operator::coulomb, libint2::BraKet::xx_xx, 0>(
                        obs[s1], obs[s2], obs[s3], obs[s4]
                    );
                    const auto* buf_1234 = buf[0];
                    for (auto f1 = 0ul, f1234 = 0ul; f1 != n1; ++f1) {
                        const auto bf1 = f1 + bf1_first;
                        for (auto f2 = 0ul; f2 != n2; ++f2) {
                            const auto bf2 = f2 + bf2_first;
                            for (auto f3 = 0ul; f3 != n3; ++f3) {
                                const auto bf3 = f3 + bf3_first;
                                for (auto f4 = 0ul; f4 != n4; ++f4, ++f1234) {
                                    const auto bf4 = f4 + bf4_first;
                                    const auto integral = buf_1234[f1234];
                                    printf("(%u, %u, %u, %u), %lf \n", bf1, bf2, bf3, bf4, integral);
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    libint2::finalize();
    return 0;
}

which I compile with

g++ -Lpath/to/libint2/lib/ -Wl,-rpath=path/to/libint2/lib/ -Ipath/to/libint2/include -Ipath/to/libint2/build/include -Ipath/to/eigen/ -o test test.cc -lint2

and link against libint2 2.7.0beta6.

When comparing libint2 integrals with results from other integral libraries such as pyscf/libcint, I've noticed that there are some mismatches in the integrals that involve some of the d functions, and in particular, some permutation symmetries are not upheld. For example, in chemists' notation:

(0, 0, 3, 8), 0.283853 // <-- matches libcint
(3, 8, 0, 0), 0.491649 // <-- incorrect by a factor of 1.732

(1, 1, 3, 6), 0.283853 // <-- matches libcint
(3, 6, 1, 1), 0.491649 // <-- incorrect by a factor of 1.732

(0, 0, 8, 8), 0.836689 // <-- matches libcint
(8, 8, 0, 0), 1.449188 // <-- incorrect by a factor of 1.732

(0, 0, 4, 4), 0.851560 // <-- matches libcint
(4, 4, 0, 0), 0.491649 // <-- incorrect by a factor of 0.577 = 1/1.732

(0, 0, 2, 2), 3.203902 // <-- matches libcint
(2, 2, 0, 0), 3.203902 // <-- no issue here

(4, 4, 8, 8), 0.602537 // <-- matches libcint
(8, 8, 4, 4), 0.602537 // <-- no issue here

From what I've understood from libint2's documentation, the order of the basis functions is

0 px
1 py
2 pz
3 dxx
4 dxy
5 dxz
6 dyy
7 dyz
8 dzz

I am just wondering if I've done something wrong when retrieving the integrals from libint2's engine, or if I have missed out an option somewhere to scale the basis functions properly. Any help is appreciated! Many thanks!

evaleev / libint

Libint2 gives some two-electron integrals that are permutation-symmetry-inconsistent? #216