Integrals

The integrals module (namespace lible::ints) provides utilities to calculate molecular integrals over Gaussian-type basis functions. The implementation is based on the McMurchie-Davidson scheme; calculation of the two-electron integrals makes use of the SHARK algorithm. Some of the notations and conventions here might bear a close resemblance to what is in the “Molecular Electronic-Structure Theory”. This book served as the main reference in implementing the library. In order to use the library, include the main header <lible/ints/ints.hpp> in your source code.

Definitions & conventions

Lible provides utilities to calculate molecular integrals over contracted Gaussian type atomic orbitals (GTAOs). A contracted GTAO can be written as

\[G^{lm}_{\mu} = N_{\mu} \sum^K_{i=1} d_{\mu,i} g^{lm}_{\mu,i}\]

where \(g^{lm}_{\mu,i}\) denotes a primitive Gaussian basis function. The contraction depth is denoted by \(K\) and the contraction coefficients by \(d_{\mu,i}\). The normalization coefficient \(N_{\mu}\) is usually obtained from the overlap integral

\[(\mu|\nu) = \int^{\infty}_{-\infty} (G^{lm}_{\mu} G^{lm}_{\nu}) d^3\mathbf{r} \; \Rightarrow \; N_{\mu} = 1.0 / \sqrt{(\mu|\mu)}\]

The function lible::ints::calcShellNorms may be used for that purpose. The primitive Gaussian basis functions are normalized, and depend on the orbital angular momentum \(l\) and its projection, \(m_l = -l, -l + 1, \ldots, l\). Omitting here the \(\mu\)-label, a primitive Gaussian basis function can be written as

\[g^{lm}_{i}(\mathbf{r}, \mathbf{a}, \mathbf{A}) = N_{l} (a_i) S_{lm} (\mathbf{r}, \mathbf{A}) e^{-a_i r^2_A}\]

The quantities in this expression are defined as:

1. \(\mathbf{r} = (x, y, z)\) – an arbitrary point in space.

2. \(\mathbf{A} = (A_x, A_y, A_z)\) – point in space where the basis function is centered, i.e., a nucleus.

3. \(\mathbf{r}_A = (x - A_x, y - A_y, z - A_z)\).

4. \(\mathbf{a} = (a_1,\ldots,a_K)\) – list of Gaussian primitive exponents.

5. \(N_{l} (a_i)\) – pure (harmonic) Gaussian primitive normalization coefficient, lible::ints::purePrimitiveNorm, that is calculated analytically from

\[N_{l}(a_i) = \sqrt{\frac{(2a_i/\pi)^{3/2}(4a_i)^l}{(2l - 1)!!}}\]

6. \(S_{lm} (\mathbf{A})\) – a real-valued solic harmonic. The explicit expressions for these shall not be given here. Essentially, the solid harmonics can be expressed in terms of Cartesian directions through the transformation \(t_{lm;ijk}\),

   \[S_{lm} = \sum_{ijk} t_{lm,ijk} x_A^i y_A^j z_A^k\]

   where the transformation coefficients are given by given by eq. (9.1.9) in “Molecular Electronic-Structure Theory”.

The Cartesian to spherical transformation is accessible from the library via the function lible::ints::sphericalTrafo which returns the transformation coefficients up to \(l=9\). Lible follows a convention where the spherical atomic orbitals are ordered alternatingly by the \(m_l\) quantum number, as in the table:

l	\(m_l\)
0	0
1	0, 1,-1
2	0, 1,-1, 2,-2
3	0, 1,-1, 2,-2, 3,-3

The number of spherical harmonic atomic orbitals is given by \(N = 2l + 1\), lible::ints::numSphericals.

The Cartesian Gaussian functions are ordered alphabetically by the Cartesian exponents. For the first few angular momenta, \(l = 0-3\), the Cartesian exponents are given by:

\(l\)	Cartesian exponents
0	1
1	(1, 0, 0), (0, 1, 0), (0, 0, 1)
2	(2, 0, 0), (1, 1, 0), (1, 0, 1), (0, 2, 0), (0, 1, 1), (0, 0, 2)
3	(3, 0, 0), (2, 1, 0), (2, 0, 1), (1, 2, 0), (1, 1, 1), (1, 0, 2), (0, 3, 0), (0, 2, 1), (0, 1, 2), (0, 0, 3)

To get the Cartesian exponents for arbitrary angular momentum, l, use lible::ints::cartExps. The total number of Cartesian exponents is given by \(N = (l + 1)(l + 2) / 2\), lible::ints::numCartesians.

Available integral kernels

The library provides kernel functions for various integral types. A kernel function calculates the integrals for one batch, which is either one shell doublet, triplet or a quartet. In the spirit of the SHARK method, all of the kernel functions return the integrals normalized and in a spherical basis. The available integral kernels are summarized in the table below:

Integral	Function
\((a|b)\)	lible::ints::overlapKernel
\(\boldsymbol{\nabla}_{AB} (a|b)\)	lible::ints::overlapD1Kernel
\((a|{-\tfrac{1}{2}}\nabla^2|b)\)	lible::ints::kineticEnergyKernel
\(\boldsymbol{\nabla}_{AB} (a|{-\tfrac{1}{2}}\nabla^2|b)\)	lible::ints::kineticEnergyD1Kernel
\((a|\sum_q \frac{-q}{r_{1q}}|b)\)	lible::ints::externalChargesKernel
\((a|\sum_q \frac{-q \cdot \text{erf}(\omega r_{1q})}{r_{1q}}|b)\)	lible::ints::externalChargesErfKernel
\(\boldsymbol{\nabla}_{AB} (a|\sum_q \frac{-q}{r_{1q}}|b)\)	lible::ints::externalChargesD1Kernel
\(\boldsymbol{\nabla}_{AB} (a|\sum_q \frac{-q \cdot \text{erf}(\omega r_{1q})}{r_{1q}}|b)\)	lible::ints::externalChargesErfD1Kernel
\(\boldsymbol{\nabla}_{q} (a|\frac{-q}{r_{1q}}|b)\)	lible::ints::externalChargesOperatorD1Kernel
\(\boldsymbol{\nabla}_{q} (a| \frac{-q \cdot \text{erf}(\omega r_{1q})}{r_{1q}}|b)\)	lible::ints::externalChargesOperatorErfD1Kernel
\((a|\frac{-q}{r_{1q}}|b)\)	lible::ints::potentialAtExternalChargesKernel
\((a| \frac{-q \cdot \text{erf}(\omega r_{1q})}{r_{1q}}|b)\)	lible::ints::potentialAtExternalChargesErfKernel
\((a|\boldsymbol{\hat{\mu}}_O|b)\)	lible::ints::dipoleMomentKernel
\(\boldsymbol{\nabla}_{A} \times \boldsymbol{\nabla}_{B} (a|\sum_q \frac{-q}{r_{1q}}|b)\)	lible::ints::spinOrbitCoupling1ElKernel
\((a|{-\boldsymbol{\nabla}}|b)\)	lible::ints::momentumKernel
\((a|{-\boldsymbol{r} \times \boldsymbol{\nabla}}|b)\)	lible::ints::angularMomentumKernel
\((a|\hat{p}_i \hat{V} \hat{p}_j|b)\)	lible::ints::pVpKernel
\((ab|\frac{1}{r_{12}}|cd)\)	lible::ints::ERI4Kernel::operator()
\((ab|\frac{1}{r_{12}}|c)\)	lible::ints::ERI3Kernel::operator()
\((a|\frac{1}{r_{12}}|b)\)	lible::ints::ERI2Kernel::operator()
\(\boldsymbol{\nabla}_{ABCD}(ab|\frac{1}{r_{12}}|cd)\)	lible::ints::ERI4D1Kernel::operator()
\(\boldsymbol{\nabla}_{ABC}(ab|\frac{1}{r_{12}}|c)\)	lible::ints::ERI3D1Kernel::operator()
\(\boldsymbol{\nabla}_{AB}(a|\frac{1}{r_{12}}|b)\)	lible::ints::ERI2D1Kernel::operator()
\((\boldsymbol{\nabla}_{AB} \otimes \boldsymbol{\nabla}_{AB})(a|\frac{1}{r_{12}}|b)\)	lible::ints::ERI2D2Kernel::operator()
\(\boldsymbol{\nabla}_{A} \times \boldsymbol{\nabla}_{B}(ab|\frac{1}{r_{12}}|cd)\)	lible::ints::ERI4SOCKernel::operator()
\(\boldsymbol{\nabla}_{A} \times \boldsymbol{\nabla}_{B}(ab|\frac{1}{r_{12}}|c)\)	lible::ints::ERI3SOCKernel::operator()

Main Interface

Typically, integral calculation in quantum chemistry programs involves choosing a molecular geometry and a basis set. From this data, shells (lible::ints::Shell) can be constructed that contain all the information required for calculating the integrals. In Lible, the shells are taken to create a special data structure called the shell pair data (lible::ints::ShellPairData). Therefore, typical integral calculation with Lible involves the shell pair data and specifying a pair of shells with an index. Graphically, the flow of data from the initial geometry and basis set to integral calculation can be illustrated as:

It should be noted here, that sometimes, the shells are transformed into the so-called shell data (lible::ints::ShellData) data structure. This data structure can be used for calculating integrals involving auxiliary basis functions, such as \((\mu\nu|P)\).

For convenience, it is possible to calculate some integrals directly, without utilizing the shell pair data. Lible provides a special data structure that records all the information about geometry, basis sets and shells for that purpose: lible::ints::Structure. Using lible::ints::Structure, integrals can be calculated such that the management of shell pair data is done under the hood. Graphically, this looks as follows:

This approach can be preferable for testing and prototyping. For large scale implementation, the previous scheme is probably preferable. For example, the function lible::ints::eri4 calculates all of the four-center two-electron integrals. For large systems, the returned 4D array of integrals would become too large to be stored in memory.

Let us assume that Lible is properly built/installed and linked against your code. To use the library for calculating integrals, include the main header #include <lible/ints/ints.hpp> in your source code. This header file constitutes the so-called main interface. The main interface contains inclusions of some other Lible header files:

In the following we shall expose the contents of the main header file. After that, the documentation of the other header files will be provided.

<lible/ints/ints.hpp>

For the code snippets shown below, assume the following data is available:

std::string basis_set = "def2-svp";
std::string basis_set_aux = "def2-universal-jkfit";
std::vector<int> atomic_nrs_h2o{8, 1, 1};
std::vector<std::array<double, 3>> coords_h2o_ang{
    {0.00000, 0.00000, 0.11779},
    {0.00000, 0.75545,-0.47116},
    {0.00000,-0.75545,-0.47116}};

vec2d overlap(const Structure &structure)

Calculates the overlap integrals. Uses OpenMP parallelization.

lible::ints::Structure structure(basis_set, atomic_nrs, coords_h2o_ang);
lible::vec2d ovlp_ints = lible::ints::overlap(structure);

vec2d overlapKernel(size_t ipair, const ShellPairData &sp_data)

Calculates a batch of overlap integrals. The shell pair index ipair refers to a pair of shells in the shell pair data sp_data.

for (size_t ipair = 0; ipair < sp_data.n_pairs_; ipair++)
{
    // The index `ipair` is used by the kernel function to read the data necessary for
    // integral calculation.
    lible::vec2d ints_batch = lible::ints::overlapKernel(ipair, sp_data);
    ...
}

std::array<vec2d, 6> overlapD1Kernel(size_t ipair, const ShellPairData &sp_data)

Calculates a batch of first derivative overlap integrals. Returns the derivative integrals for the atomic centers involved in the shell pair (ipair): \(({\partial_A}_x, {\partial_A}_y, {\partial_A}_z, {\partial_B}_x, {\partial_B}_y, {\partial_B}_z)\).

for (size_t ipair = 0; ipair < sp_data.n_pairs_; ipair++)
{
    std::array<lible::vec2d, 6> ints_batch = lible::ints::overlapD1Kernel(ipair, sp_data);

    // The indices of the atoms (A and B) involved in the integrals batch can be obtained
    // from the shell pair data.
    size_t iatom_a = atomic_idxs_[2 * ipair];
    size_t iatom_b = atomic_idxs_[2 * ipair + 1];
}

vec2d kineticEnergy(const Structure &structure)

Calculates the kinetic energy integrals. Uses OpenMP parallelization.

vec2d kineticEnergyKernel(size_t ipair, const ShellPairData &sp_data)

Calculates a batch of kinetic energy integrals.

std::array<vec2d, 6> kineticEnergyD1Kernel(size_t ipair, const ShellPairData &sp_data)

Calculates a batch of first derivative kinetic energy integrals.

vec2d nuclearAttraction(const Structure &structure)

Calculates nuclear attraction integrals. Uses OpenMP parallelization.

vec2d nuclearAttractionErf(const Structure &structure, const std::vector<double> &omegas)

Calculates the attenuated Coulomb attraction integrals. Uses OpenMP parallelization.

Important

The number of attenuation parameters omegas must equal the number of atoms in the structure.

vec2d externalCharges(const std::vector<std::array<double, 4>> &point_charges, const Structure &structure)

Calculates the one-electron Coulomb integrals with given point charges. Uses OpenMP parallelization.

vec2d externalChargesErf(const std::vector<std::array<double, 4>> &point_charges, const std::vector<double> &omegas, const Structure &structure)

Calculates the attenuated one-electron Coulomb integrals with given point charges. Uses OpenMP parallelization.

Important

The numbers of attenuation parameters omegas and point charges point_charges must be equal.

vec2d externalChargesKernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of one-electron Coulomb integrals for a list of given charges.

for (const lible::ints::ShellPairData &sp_data : sp_data_all)
{
    // The boys grid must be initialized with the angular momentum correct angular momentum:
    auto [la, lb] = sp_data.getLPair();
    int lab = la + lb;
    lible::ints::BoysGrid boys_grid(lab);

    for (size_t ipair = 0; ipair < sp_data.n_pairs_; ipair++)
    {
        lible::vec2d ints_batch = lible::ints::externalChargesKernel(
            ipair, charges, boys_grid, sp_data);

        // Consume the integrals..
    }
}

vec2d externalChargesErfKernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const std::vector<double> &omegas, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of attenuated Coulomb integrals.

Important

The numbers of attenuation parameters omegas and charges charges must be equal.

std::array<vec2d, 6> externalChargesD1Kernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of first-derivative one-electron Coulomb integrals.

Important

Due to differentiation, the Boys function grid has to be initialized with total angular momentum incremented by one:

auto [la, lb] = sp_data.getLPair();
int lab1 = la + lb + 1;
lible::ints::BoysGrid boys_grid(lab1);

std::array<vec2d, 6> externalChargesErfD1Kernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const std::vector<double> &omegas, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of first-derivative attenuated one-electron Coulomb integrals.

Important

Due to differentiation, the Boys function grid has to be initialized with total angular momentum incremented by one:

auto [la, lb] = sp_data.getLPair();
int lab1 = la + lb + 1;
lible::ints::BoysGrid boys_grid(lab1);

Important

The number of the attenuation parameters omegas must equal the number of charges, \((x, y, z, q)\), charges.

std::vector<std::array<vec2d, 3>> externalChargesOperatorD1Kernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of one-electron Coulomb operator derivative integrals. The returned list of integrals has the length of the given charges, charges.

std::vector<std::array<vec2d, 3>> externalChargesOperatorErfD1Kernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const std::vector<double> &omegas, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of attenuated one-electron Coulomb operator derivative integrals. The returned list of integrals has the length of the given charges, charges.

Important

The Boys function grid has to be initialized with total angular momentum of \(l = l_a + l_b + 1\).

Important

The number of the attenuation parameters omegas must equal the number of charges, \((x, y, z, q)\), charges.

std::vector<vec2d> potentialAtExternalChargesKernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of one-electron Coulomb integrals at the given charges, \((x, y, z, q)\). Returns a list of integrals with the length of given charges.

std::vector<vec2d> potentialAtExternalChargesErfKernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const std::vector<double> &omegas, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of attenuated one-electron Coulomb integrals at the given charges, \((x, y, z, q)\). Returns a list of integrals with the length of given charges.

Important

The list of charges must have the same length as the attenuated parameters (omegas).

std::array<vec2d, 3> dipoleMoment(const std::array<double, 3> &origin, const Structure &structure)

Calculates dipole moment integrals. Uses OpenMP parallelization.

std::array<vec2d, 3> dipoleMomentKernel(size_t ipair, const std::array<double, 3> &origin, const ShellPairData &sp_data)

Calculates a batch of dipole moment integrals.

std::array<vec2d, 3> spinOrbitCoupling1El(const Structure &structure)

Calculates spin-orbit coupling (SOC) one-electron integrals. Uses OpenMP parallelization.

std::array<vec2d, 3> spinOrbitCoupling1ElKernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of of spin-orbit coupling (SOC) one-electron integrals.

Important

The boys grid must be initialized with \(l = l_a + l_b + 1\).

std::array<vec2d, 3> momentum(const Structure &structure)

Calculates linear momentum integrals. Uses OpenMP parallelization.

std::array<vec2d, 3> momentumKernel(size_t ipair, const ShellPairData &sp_data)

Calculates a batch of linear momentum integrals.

std::array<vec2d, 3> angularMomentum(const std::array<double, 3> &origin, const Structure &structure)

Calculates angular momentum integrals. Uses OpenMP parallelization.

std::array<vec2d, 3> angularMomentumKernel(size_t ipair, const std::array<double, 3> &origin, const ShellPairData &sp_data)

Calculates a batch of angular momentum integrals.

arr2d<vec2d, 3, 3> pVpIntegrals(const Structure &structure)

Calculates the \(\hat{p}_i \hat{V}\hat{p}_j\) integrals. Used in X2C: https://doi.org/10.1063/1.4803693. Uses OpenMP parallelization.

arr2d<vec2d, 3, 3> pVpKernel(size_t ipair, const std::vector<std::array<double, 4>> &charges, const BoysGrid &boys_grid, const ShellPairData &sp_data)

Calculates a batch of \(\hat{p}_i \hat{V}\hat{p}_j\) integrals.

Important

The boys grid must be initialized with \(l = l_a + l_b + 2\).

std::vector<double> eri2Diagonal(const Structure &structure)

Calculates the diagonal part of the two-center Coulomb repulsion integrals in auxiliary basis: \((K|K)\). Uses OpenMP parallelization.

Important

The structure must be initialized with an auxiliary basis set beforehand.

vec2d eri4Diagonal(const Structure &structure)

Calculates the diagonal part of the four-center Coulomb repulsion integrals, \((\mu\nu|\mu\nu)\). Uses OpenMP parallelization.

vec3d eri3(const Structure &structure)

Calculates the three-center Coulomb repulsion integrals, \((\mu\nu|K)\). Uses OpenMP parallelization.

vec4d eri4(const Structure &structure)

Calculates the four-center Coulomb repulsion integrals, \((\mu\nu|\kappa\tau)\). Uses OpenMP parallelization

BasisAtom basisForAtom(int atomic_nr, const std::string &basis_set)

Returns the main basis set for for an atom.

BasisAtom basisForAtomAux(int atomic_nr, const std::string &aux_basis_set)

Returns the auxiliary basis set for an atom.

basis_atoms_t basisForAtoms(const std::vector<int> &atomic_nrs, const std::string &basis_set)

Returns the main basis set for the given atoms.

basis_atoms_t basisForAtomsAux(const std::vector<int> &atomic_nrs, const std::string &aux_basis_set)

Returns the auxiliary basis set for the given atoms.

std::vector<std::tuple<int, int, double>> sphericalTrafo(int l)

Returns the Cartesian to spherical transformation for the given angular momentum.

std::vector<double> atomic_orbitals_cart(lible::ints::numCartesians(l));

// Calculate values of atomic orbitals

auto spherical_trafo = lible::ints::sphericalTrafo(l);
std::vector<double> atomic_orbitals_sph(lible::ints::numSphericals(l));
// Transform AOs in Cartesian basis to spherical basis.
for (const auto& [mu, mu_, val] : spherical_trafo)
    atomic_orbitals_sph[mu] += val * atomic_orbitals_cart[mu_];

<lible/ints/shell.hpp>

struct BasisShell

Structure representing the basis set for an atomic orbital shell.

int l_

Angular momentum of the shell.

std::vector<double> exps_

Gaussian primitive exponents.

std::vector<double> coeffs_

Contraction coefficients of the Gaussian primitives.

using basis_shells_t = std::vector<lible::ints::BasisShell>

Type for representing basis sets on a list of shells.

struct BasisAtom

Structure representing the basis set on an atom.

int atomic_nr_

Atomic number of the atom.

basis_shells_t basis_shells_

List of shell basis sets on the atom.

using basis_atoms_t = std::vector<BasisAtom>

Type representing a list of basis sets on atoms. Object of this type can be used to represent the entire molecular basis set (main or auxiliary).

struct Shell

Structure representing an atomic orbital shell. Contains essential data for calculating integrals. rambleramble.

int l_

Angular momentum of the shell.

int z_

Atomic number of the shell’s atom.

size_t dim_cart_

Number of Cartesian atomic orbitals.

size_t dim_sph_

Number of spherical atomic orbitals.

size_t ofs_cart_

Starting position of the current shell’s Cartesian atomic orbitals in the list of all AOs.

size_t ofs_sph_

Starting position of the current shell’s spherical atomic orbitals in the list of all AOs.

size_t idx_

Index of the current shell in the list of all shells.

size_t idx_atom_

Index of the current shell’s atom in the list of all atoms.

std::array<double, 3> xyz_coords_

Cartesian coordinates of the shell’s atom.

std::vector<double> exps_

Gaussian primitive exponents.

std::vector<double> coeffs_

Contraction coefficients of the Gaussian primitives.

std::vector<double> norms_

Normalization coefficients of the spherical atomic orbitals.

std::vector<double> norms_prim_

Normalization coefficients of the Gaussian primitives.

<lible/ints/shell_pair_data.hpp>

struct ShellData

struct ShellPairData