Skip to content

using modern CPlusPlus API

Jump to bottom
Eduard Valeyev edited this page Dec 11, 2023 · 21 revisions

intro

This document describes the features of the "modern" C++ libint interface (version 2.6.0). To get a quick feel the interface see:

  1. a simple direct closed-shell Hartree-Fock energy code: tests/hartree-fock/hartree-fock.cc
  2. a more elaborate (shared-memory) parallel code for closed-shell Hartree-Fock energy (+ forces and core hessian (i.e. sans orbital response)) with efficiency comparable to production codes: tests/hartree-fock/hartree-fock++.cc

prerequisites

To be able to use the C++11 interface you will need to have the following:

  1. C++ compiler that supports C++11 standard. C++11 standard is the third most recent C++ international standard for C++, hence most modern compilers support it fully. Most compilers will default to C++11 or later, but if yours does not you will need to use -std=c++11 compiler flag to turn this support on.
  2. libint library, configured (cmake ...) and installed (cmake --target install...). To get started you could try this library . In cmake you can just do find_package(Libint2); if compiling a program that depends on libint manually you will need to specify the location of libint headers (-I/libint_prefix/include, where libint_prefix is the installation prefix specified when you configured libint; by default it's /usr/local). You will also need to specify the location of libint library (-L/libint_prefix/lib).
  3. Eigen C++ library: a headers-only (i.e. no need to compile) linear algebra library. Download, unpack, and pass the location to the compiler (-I/path/to/eigen/directory).

quick tour of API

to compute integrals you need to

  1. include libint2.hpp header
  2. bootstrap libint (this step initializes some internal data)
  3. create a basis set object (potentially, more than one if you are doing, for example, density-fitting)
  4. create an integral engine
  5. iterate over shell tuples (shell sets), compute integrals, and digest them (store, contract, print, etc.)

boostrap libint

To bootstrap libint call initialize(). When you are done, call finalize(). Both of these functions (and other libint-specific classes) live in namespace libint2:

#include <libint2.hpp>

int main(int argc, char* argv[]) {
  libint2::initialize();  // safe to use libint now .. do `libint2::initialize(true)` to produce diagnostic messages

  // all other code snippets go here

  libint2::finalize();  // do not use libint after this

  // can repeat the libint2::initialize() ... finalize() cycle as many times as
  // necessary

  return 0;
}

For brevity I will omit libint2:: and std::: this assumes you added using namespace libint2; and using namespace std; (NB this is not a good idea to add these statements in global scope of C++ header files).

create a basis set

Here's how to create a basis set object using the built-in basis set library

// #include <fstream>, <string>, and <vector>

string xyzfilename = "/path/to/input_dot_xyz_file"; // see XYZ format description at http://en.wikipedia.org/wiki/XYZ_file_format
ifstream input_file(xyzfilename);
vector<Atom> atoms = read_dotxyz(input_file);

BasisSet obs("cc-pVDZ", atoms);

Think of BasisSet as a slightly decorated std::vector<Shell>. This means that you can manipulate it as a std::vector directly, e.g. append shells, reorder shells, sort shells, etc. For example, here's how to sort shells by their angular momenta:

// #include <algorithm> and <iterator>

// print out the original basis
std::copy(begin(obs), end(obs),
          std::ostream_iterator<Shell>(std::cout, "\n"));

// stable sort preserves the original order where possible
stable_sort(begin(obs), end(obs),
            [](const Shell& a, const Shell& b){
               return a.contr[0].l < b.contr[0].l;
             });

// print out the resorted basis
std::copy(begin(obs), end(obs),
          std::ostream_iterator<Shell>(std::cout, "\n"));

NB this example assumes a basis set with segmented or simple general contraction scheme, where each shell has contractions of 1 angular momentum only. Note that the SP shells in Pople-style basis sets are split to make a segmented basis.

BasisSet exists to make basis set construction easier (e.g. by reading standard basis sets from a basis set stash accompanying libint). However C++11's uniform initialization syntax makes it very easy to construct shells and basis sets directly, as shown in tests/hartree-fock/hartree-fock.cc .

BasisSet assumes that all shells of angular momentum > 1 are solid harmonics Gaussians, with few exceptions as denoted in Gaussian09 manual. It is trivial to adjust the type of shells after constructing the basis; the following example shows how to force the use cc-pVTZ basis with Cartesian shells for all d functions:

BasisSet obs("cc-pVTZ", atoms);
for(auto& shell: obs)  // NB range-based for
  for(auto& contraction: shell.contr)
    if (contraction.l == 2) contraction.pure = false;

or even simpler:

BasisSet obs("cc-pVTZ", atoms);
obs.set_pure(false); // use cartesian gaussians

create an integral engine

Prior to version 2.2 libint provided separate engine types for 1-body (OneBodyEngine) and 2-body (TwoBodyEngine) operators. No longer: a single engine type (Engine) handles all integrals. Here's how to create an integral engine for computing overlap integrals:

Engine s_engine(Operator::overlap,  // will compute overlap ints
                obs.max_nprim(),    // max # of primitives in shells this engine will accept
                obs.max_l()         // max angular momentum of shells this engine will accept
               );

The first parameter to the constructor of Engine is a value of enum class Operator type that specifies the operator (actually, a set of one or more operators; see below) that the engine will evaluate. The currently supported values are listed in the following table.

Operator set Description Parameters (= default)
Operator::overlap overlap none
Operator::kinetic kinetic energy, $-\frac{1}{2} \nabla^2$ none
Operator::nuclear (negative) Coulomb potential due to point charges, $- \sum\limits_{I=1}^{\text{# of charges}} Z_I / |\vec{r} - \vec{R}_I|$ std::vector<std::pair<Real, std::array<Real, 3>>> = {}, vector of ${Z_I,\vec{R}_I}$ pairs
Operator::emultipole1 overlap + 3 components, {x, y, z}, of the (Cartesian) electric dipole operator std::array<Real, 3> = {0, 0, 0}, multipole origin
Operator::emultipole2 operators in Operator::emultipole1 + 6 components, {xx, xy, xz, yy, yz, zz}, of the Cartesian electric quadrupole operator same as Operator::emultipole1
Operator::emultipole3 operators in emultipole2 + 10 components, {xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz}, of the Cartesian electric octopole operator same as Operator::emultipole1
Operator::opVop small-component nuclear potential, $\left(\sigma \cdot \hat{\bf p} \right) V \left(\sigma \cdot \hat{\bf p} \right)$, in Pauli representation as 4 components, { $\hat{\bf p} \cdot V \hat{\bf p}$, $(\hat{\bf p} \times V \hat{\bf p})_z$, $(\hat{\bf p} \times V \hat{\bf p})_x$, $(\hat{\bf p} \times V \hat{\bf p})_y$} same as Operator::nuclear
Operator::sphemultipole real spherical electric multipole moments, $O_{l,m} \equiv \mathcal{N}^{\text{sign}(m)}_{l,|m|}$ where $\mathcal{N}^{\pm}_{l,m}$ is defined in J.M. Pérez-Jordá and W. Yang, J Chem Phys 104, 8003 (1996) . To obtain the traditional real solid harmonics $C^m_l$ and $S^m_l$ defined here multiply these harmonics by $(-1)^m \sqrt{(2 - \delta_{m,0}) (l + |m|)! (l - |m|)!}$ . same as Operator::emultipole1
Operator::coulomb two-body Coulomb, $r_{12}^{-1} \equiv |\vec{r}_1 - \vec{r}_2|^{-1}$ none
Operator::cgtg contracted Gaussian geminal, $f^{\text{cGTG}}{12} \equiv \sum\limits{i=1}^{n} c_i \exp(-\alpha_i r_{12}^2)$ std::vector<std::pair<Real, Real>> = {}, vector of ${c_i, \alpha_i}$ pairs
Operator::cgtg_x_coulomb contracted Gaussian geminal times Coulomb, $f^{\text{cGTG}}{12}/r{12}$ same as Operator::cgtg
Operator::delcgtg2 gradient norm of contracted Gaussian geminal, $|\nabla \dot f^{\text{cGTG}}_{12}|^2$ same as Operator::cgtg
Operator::stg Slater geminal, $f^{\text{STG}}{12} \equiv \exp(-\zeta r{12})$ Real = 0, geminal exponent parameter $\zeta$
Operator::yukawa Slater geminal times Coulomb, $f^{\text{STG}}{12}/r{12}$ same as Operator::stg
Operator::r12 the anti-Coulomb kernel, $r_{12} \equiv |\vec{r}_1 - \vec{r}_2|$ none
Operator::erf_coulomb the erf-attenuated Coulomb kernel, $\text{erf}(\omega r_{12})/r_{12}$ Real = 0, attenuation parameter $\omega$
Operator::erfc_coulomb the erfc-attenuated Coulomb kernel, $\text{erfc}(\omega r_{12})/r_{12}$ Real = 0, attenuation parameter $\omega$
Operator::erf_nuclear the erf-attenuated Coulomb potential due to point charges, $- \sum\limits_{I=1}^{\text{# of charges}} Z_I \text{erf}\left( |\vec{r} - \vec{R}_I| \right)/ |\vec{r} - \vec{R}_I|$ std::tuple<Real, std::vector<std::pair<Real, std::array<Real, 3>>>> = {0,{}}, attenuation parameter + vector of ${Z_I,\vec{R}_I}$ pairs
Operator::erfc_nuclear the erfc-attenuated Coulomb potential due to point charges, $- \sum\limits_{I=1}^{\text{# of charges}} Z_I \text{erfc}\left( |\vec{r} - \vec{R}_I| \right)/ |\vec{r} - \vec{R}_I|$ same as Operator::erf_nuclear

Note that if you want to compute integrals over multiple basis sets, the second and third params (max_prim and max_l) must be the maximum values of all basis sets that you need to use. These parameters will determine how much memory each engine will reserve.

In this example only 3 parameters are provided to the constructor; these parameters do not have defaults and can always be specified. Several more parameters can be set to non-default values. The most important is the fourth parameter that specifies the order of geometric derivatives. Its default value is 0, i.e. the engine will compute non-differentiated integrals by default. The maximum allowed value is 2 (can be increased, if needed: please contact the developer). Other defaulted parameters will be illustrated below.

Making engines for some types of integrals may require additional parameters to define the operator. For example, nuclear attraction integrals (integrals of the Coulomb potential due to the point charges) need to know where the charges are located. This information can be provided to the constructor, or by using Engine::set_params() method:

// this engine will compute nuclear attraction ints
// by default there are no charges
Engine v_engine(Operator::nuclear,
                obs.max_nprim(), obs.max_l());
v_engine.set_params(make_point_charges(atoms));  // convert `atoms` to point charges
                                                 // classical charges in QM/MM need extra work

Another example are electric (Cartesian) multipole integrals. The default is to use the coordinate system origin, {0,0,0}, to define the multipole expansion. If a different point is to be used for the multipole expansion its Cartesian coordinates must be provided, as an array<double,3>, via set_params:

Engine mu_engine(Operator::emultipole1, // will compute overlap + electric dipole moments
                 obs.max_nprim(), obs.max_l());
mu_engine.set_params(array<double,3>{2.0, 3.0, 4.0}); // use x=2,y=3,z=4 as the multipole origin

Making engines for two-body integrals is equally simple:

Engine eri_engine(Operator::coulomb, obs.max_nprim(), obs.max_l());

Several other two-body kernel types are supported. For example, range-separated DFT methods require integrals over erf(c)-attenuated Coulomb kernel:

double omega = 1.2;
Engine erfcoulomb_engine(Operator::erf_coulomb, obs.max_nprim(), obs.max_l());
erfcoulomb_engine.set_params(omega);

Exponentially-attenuated Coulomb (aka Yukawa) kernel is also supported as of version 2.6.0:

double gamma = 1.0;
Engine yukawa_engine(Operator::yukawa, obs.max_nprim(), obs.max_l());
yukawa_engine.set_params(gamma);

Explicitly correlated (R12/F12) methods typically use exponential (Slater-type) correlation kernel, whose integrals are available in libint 2.6.0 and later:

double gamma = 1.0;
Engine stg_engine(Operator::stg, obs.max_nprim(), obs.max_l());
stg_engine.set_params(gamma);

Prior to 2.6.0 integrals over Slater-type geminal kernel were computed by expanding it into Gaussian geminals. Several integrals over general contracted Gaussian geminal are supported and can be used to evaluate integrals over suitably-well-behaved spherically-symmetric geminals. For example, to compute integrals over a contracted Gaussian-Type Geminals (cGTG), f(r12) = 0.1 exp(-0.2 r122) - 0.3 exp(-0.4 r12) + 0.5 exp(-0.6 r12), use this code:

vector<pair<double,double>> cgtg_params{{0.2,0.1}, {0.4,-0.3}, {0.6, 0.5}};
Engine cgtg_engine(Operator::cgtg, obs.max_nprim(), obs.max_l());
cgtg_engine.set_params(cgtg_params);
// an alternative is to provide operator params to ctor directly
// but all ctor arguments preceeding params must be given explicitly:
// Engine cgtg_engine(Operator::cgtg, obs.max_nprim(), obs.max_l(),
//                    0,   // regular ints, i.e. not derivs
//                    numeric_limits<double>::epsilon(),  // precision
//                    cgtg_params);

Other kernel types derived from cGTG are used similarly:

Engine cgtg_coulomb_engine(obs.max_nprim(), obs.max_l());
cgtg_coulomb_engine.set_params(cgtg_params);

compute integrals

Using the engines to compute integrals is trivial. Here's how to compute the entire set of overlap integrals:

auto shell2bf = obs.shell2bf(); // maps shell index to basis function index
                                // shell2bf[0] = index of the first basis function in shell 0
                                // shell2bf[1] = index of the first basis function in shell 1
                                // ...
const auto& buf_vec = s_engine.results(); // will point to computed shell sets
                                          // const auto& is very important!

for(auto s1=0; s1!=obs.size(); ++s1) {
  for(auto s2=0; s2!=obs.size(); ++s2) {

    cout << "compute shell set {" << s1 << "," << s2 << "} ... ";
    s_engine.compute(obs[s1], obs[s2]);
    cout << "done" << endl;
    auto ints_shellset = buf_vec[0];  // location of the computed integrals
    if (ints_shellset == nullptr)
      continue;  // nullptr returned if the entire shell-set was screened out

    auto bf1 = shell2bf[s1];  // first basis function in first shell
    auto n1 = obs[s1].size(); // number of basis functions in first shell
    auto bf2 = shell2bf[s2];  // first basis function in second shell
    auto n2 = obs[s2].size(); // number of basis functions in second shell

    // integrals are packed into ints_shellset in row-major (C) form
    // this iterates over integrals in this order
    for(auto f1=0; f1!=n1; ++f1)
      for(auto f2=0; f2!=n2; ++f2)
        cout << "  " << bf1+f1 << " " << bf2+f2 << " " << ints_shellset[f1*n2+f2] << endl;
  }
}

Each call to Engine::compute returns a shell-block of integrals, packed into a linear array. Computing one-shell-block-at-a-time is done for the sake of efficiency.

What about dipole integrals, for example? Calling compute method on mu_engine will compute multiple (namely, four) shell blocks, one block per corresponding operator in the integral set computed by that engine. The 4 shell blocks correspond to overlap followed by the integrals over x, y, and z compoments of the electric dipole operator. Here's a code snippet that shows how to digest the integrals:

const auto& buf_vec = mu_engine.results(); // will point to computed shell sets
                                           // const auto& is very important!
for(auto s1=0; s1!=obs.size(); ++s1) {
  for(auto s2=0; s2!=obs.size(); ++s2) {

    cout << "compute shell set {" << s1 << "," << s2 << "} ... ";
    mu_engine.compute(obs[s1], obs[s2]);
    cout << "done" << endl;

    auto bf1 = shell2bf[s1];  // first basis function in first shell
    auto n1 = obs[s1].size(); // number of basis functions in first shell
    auto bf2 = shell2bf[s2];  // first basis function in second shell
    auto n2 = obs[s2].size(); // number of basis functions in second shell

    auto s_shellset = buf_vec[0];  // skipped nullptr check for simplicity
    auto mu_x_shellset = buf_vec[1];
    auto mu_y_shellset = buf_vec[2];
    auto mu_z_shellset = buf_vec[3];
    cout << "e.g. z dipole integrals:" << endl;
    for(auto f1=0; f1!=n1; ++f1)
      for(auto f2=0; f2!=n2; ++f2)
        cout << "  " << bf1+f1 << " " << bf2+f2 << " " << mu_z_shellset[f1*n2+f2] << endl;
  }
}

Computing shell-sets of multiple related operators together is again done for the sake of efficiency.

Similarly, an engine constructed with Operator::emultipole2 will return 10 shell blocks at a time. The first 4 are the same as in the emultipole1 case, followed by the 6 unique components of the Cartesian quadrupole operators: xx, xy, xz, yy, yz, and zz.

It is sometimes conventient to evaluate target integrals multiplied by some constant scalar factor. It is possible to pass the scaling factor to the engine via Engine::prescale_by() before calling `compute():

mu_engine.prescale_by(2);
assert(mu_engine.prescaled_by() == 2);
// now mu_engine.compute() returns target integrals scaled by 2

The engine's scaling factor can be reset as often as necessary. Copies on an engine inherit its scaling factor. Note that scaling of the integrals happens before screening (hence prescale), so it may be needed to adjust the precision accordingly.

NB The method by which the integrals are returned has changed in version 2.2 in order to minimize the number of copies performed by Engine::compute internally. Instead of returning a pointer to the beginning of a memory region containing the consecutive shell-sets of integrals, compute sets the pointers to each shell-set in the vector of pointers the reference to which is obtained by Engine::results(). If all integrals are zero, the first element of this vector will be nullptr (again, this is more efficient than filling in memory block with zeroes).

API in more detail

Normalization and layout conventions

By default, libint follows a revised Common Component Architecture (CCA) standard for normalization and layout conventions; the base version of the CCA standard is described in Appendix B of J. Comp. Chem. 29, 562 (2008). Although libint libraries can be configured to use non-default orderings for Cartesian and solid harmonic Gaussians, for example to be used in applications that do not use the CCA standard, the following description with focus primarily on the CCA standard.

Cartesian Gaussian ordering and normalization

  • In the CCA standard a Cartesian shell contains basis functions in lexicographic ordering. The ordering for the four smallest shells is illustrated in the table below.
L Label Layout
0 s s
1 p x, y, z
2 d x2, xy, xz, y2, yz, z2
3 f x3, x2y, x2z, xy2, xyz, xz2, y3, y2z, yz2, z3
  • Nonstandard (non-CCA) orderings of Cartesian shells can be specified by giving configure the --with-cartgauss-ordering option when configuring the libint compiler, before generating a libint library;
  • The contraction coefficients of published basis sets usually are set so that at least some basis functions are normalized to unity. However, since some basis sets are not normalized by default, the default behavior is to rescale contraction coefficients to ensure unit normalization of all solid harmonics Gaussians and the axis-aligned (xl, yl, zl) Cartesian Gaussians. To change the default behavior call Shell::do_enforce_unit_normalization(false) before using libint.
  • The default (CCA) normalization convention for Cartesian Gaussians specifies the same normalization constant for every Gaussian in a Cartesian shell. This means that, for examle, in a d shell norms of functions {x2, y2, z2} differ from those of the {xy, xz, yz} functions. Specifically, the ${n_x, n_y, n_z }$ component of a Cartesian shell needs to be multiplied by $\sqrt{\frac{(2(n_x + n_y + n_z)-1)!!}{(2n_x - 1)!!(2n_y - 1)!!(2n_z - 1)!!}}$. This means that by default only the axis-aligned components of a Cartesian shell (i.e. with all quanta aligned on one axis, e.g. x3, y3, and z3 in an f shell) will be normalized to unity. To force every Cartesian Gaussian in a shell to have same norm call Engine::set(CartesianShellNormalization::uniform) before calling compute.

Solid Harmonic Gaussians ordering and normalization

  • The definition of unit-normalized real solid harmonics $R^{\pm}_{l,|m|}$ used by libint was taken from Int. J. Quantum Chem. 54, 83 (1995); it coincides with the definitions of the real solid harmonics elsewhere, e.g. in [WikiPedia]https://en.wikipedia.org/wiki/Solid_harmonics#Real_form). N.B. The complex solid harmonics in the IJQC paper do not include the traditional (for the quantum mechanics) Condon-Shortley phase (compare this list to Table I in the IJQC paper).
  • The CCA standard orders the solid harmonics with $R^{-}_{l,|m|}$ appearing first with |m| decreasing from l to 1, followed by $R^{+}_{l,|m|}$ with |m| increasing from 0 to l: $R^-{l,l}, R^-{l,l-1} \dots R^-{l,1}, R^+{l,0} R^+{l,1} \dots R^+{l,l}$. For example, the following table lists the CCA-ordered solid harmonic shells as linear combinations of Cartesian Gaussians, with Cartesians normalized according to the default axis-aligned convention described above:
L CCA Solid Harmonic Shell
0 1
1 $y, z, x$
2 $\sqrt{3}xy, \sqrt{3}yz, z^2 - (x^2+y^2)/2, \sqrt{3}xz, (x^2-y^2)/\sqrt{2}$
3 $ \sqrt{10} (3 x^2 y - y^3)/4, \sqrt{15} xyz, \sqrt{6} ( yz^2 - (x^2 y + y^3)/4 ), z^3 - 3 (x^2 z + y^2 z)/2, \sqrt{6} ( xz^2 - (x^3 + xy^2)/4 ), \sqrt{15} (x^2z - y^2z)/2, \sqrt{10} ( x^3 - 3 x y^2)/4$
  • Nonstandard (non-CCA) orderings of solid harmonics shells can be specified by giving configure the --with-shgauss-ordering option when configuring the libint compiler, before generating a libint library. The only non-CCA choice right now is the ordering used in the Gaussian program, defined as follows: $R^{+}{l,0}, R^{+}{l,1}, R^{-}{l,1}, R^{+}{l,2}, R^{-}_{l,2}, \dots$
  • Every function in a solid harmonic shell is normalized to unity.

Shell-set conventions:

Shell-sets (or, shell blocks) are packed similarly to how matrices are packed in row-major order, i.e. the last shell set is least significant. Thus a (p|p) shell-set of overlap integrals has the following layout: (x|x), (x|y), (x|z), (y|x), (y|y), (y|z), (z|x), (z|y), (z|z). Integrals use chemists (Mulliken) convention for function ordering. E.g. for standard 2-body integrals function 1 and 2 depend on the position of particle 1 and belong to Dirac bra and ket, respectively; functions 3 and 4 similarly depend on particle 2 position, etc.

Operator conventions:

  • Sizes and layouts of operator sets are specific to each operator set, and are currently known at compile time (this may change in the future). The layout of operator sets currently supported by libint is shown in the table below.
Operator set Layout
emultipole1 overlap, μx, μy, μz
emultipole2 overlap, μx, μy, μz, qxx, qx, qxz, qyy, qyz, qzz
emultipole3 overlap, μx, μy, μz, qxx, qx, qxz, qyy, qyz, qzz, oxxx, oxxxy, oxxz, oxyy, oxyz, oxzz, oyyy, oyyz, oyzz, ozzz
sphemultipole This operator set includes spherical electric multipole moments, $\mathcal{N}^{\pm}_{l,m}$, of order up to $l_{\rm max}$ (for a total of $(l_{\rm max}+1)^2$ operators), in the order of increasing l, with the operators of same l but different m ordered according to the solid harmonics Gaussian ordering described above. For example, for the CCA standard solid harmonics ordering the operators will appear in the following order: $\mathcal{N}^+{0,0} , \mathcal{N}^-{1,1}, \mathcal{N}^+{1,0}, \mathcal{N}^+{1,1}, \mathcal{N}^-{2,2}, \mathcal{N}^-{2,1}, \mathcal{N}^+{2,0}, \mathcal{N}^+{2,1}, \mathcal{N}^+_{2,2} \dots$
opVop p · V p, (p ⤬ V p)z, (p ⤬ V p)x, (p ⤬ V p)y
  • Shell-sets of operator sets laid out so that the operator index is most significant. E.g. the (p|p) shell-set of operator set emultipole1 is laid out with all overlap integrals first (shell-set 0), then all x-dipole integrals (shell-set 1), etc.: (x|x), (x|y), (x|z), (y|x), (y|y), (y|z), (z|x), (z|y), (z|z), (x|μx|x), (x|μx|y), (x|μx|z), (y|μx|x), (y|μx|y), (y|μx|z), (z|μx|x), (z|μx|y), (z|μx|z), (x|μy|x), (x|μy|y), (x|μy|z), (y|μy|x), (y|μy|y), (y|μy|z), (z|μy|x), (z|μy|y), (z|μy|z), (x|μz|x), (x|μz|y), (x|μz|z), (y|μz|x), (y|μz|y), (y|μz|z), (z|μz|x), (z|μz|y), (z|μz|z).

Geometric derivatives

  • libint currently supports computation of arbitrary-order geometric derivative integrals; this interface has supported arbitrary-order derivatives since version 2.7.0 (versions 2.6.0 and older only exposed up to 2nd-order derivatives of 2-body integrals).
  • N.B. geometric derivatives denote derivatives of the integral with respect to the shell origins; total derivatives with respect to geometric coordinates in the case where multiple shells share the same center can be computed trivially.
  • The number of derivatives depends on the derivative order, number of Gaussian functions in the integral, and whether the operator is depends on geometric coordinates or not:
    • Position-independent operators (all two-body operators, plus overlap and kinetic energy operators) and well as some position-dependent operators (multipole operators): the number of derivatives of order n of an integral over k Gaussians is computed as (3k) * (3k+1) * ... * (3*k+n-1) / n!. For example:
      • number of 1st derivatives of 2-center integrals = 6, ordered as: d/dx1 d/dy1, d/dz1, d/dx2, d/dy2, d/dz2, where d/dx1 is the derivative with respect to the x coordinate of the first shell in the integral
      • number of 2nd derivatives of 2-center integrals = 21, ordered as: d2/dx1dx1, d2/dx1dy1, d2/dx1dz1, d2/dx1dx2, d2/dx1dy2, d2/dx1dz2, d2/dy1dy1, d2/dy1dz1, d2/dy1dx2, ... (see Appendix B of J. Comp. Chem. 29, 562 (2008) for more details).
      • number of 2nd derivatives of 4-center integrals = 78
    • For some position-dependent operators (1-body Coulomb operator): geometric derivaties with respect to the operator coordinates are also necessary, hence the number of derivatives of order n of an integral over k Gaussians and operator dependent on l points is computed as (3*(k+l)) * (3*(k+l)+1) * ... * (3*(k+l)+n-1) / n!. For example:
      • number of 1st derivatives of nuclear attraction integrals in a molecule with 10 atoms = 36.
      • number of 2nd derivatives of nuclear attraction integrals in a molecule with 10 atoms = 666. Note that many of these derivatives are zero due to the fact that each term in the 1-body Coulomb operator depends on position of a single center only, hence the mixed derivatives with respect to coordinates of different atoms are in general zero. However, the general layout is used to support integrals over operators that depend on pairs and higher-order tuples of centers.
  • Derivative index is most significant, hence the first derivatives of operator set emultipole1 are ordered as follows: shell-set 0 contains the derivative with respect to the x coordinate of the first shell (d/dx1) of the overlap integrals, shell-set 1 -- d/dx1 of the μx integrals, ..., shell-set 4 -- d/dy1 of the overlap integrals, shell-set 23 -- d/dz2 of the μz integrals.

2- and 3-center 2-body integrals

The standard 2-body integrals include 1 functions / particle in bra and ket. The 3-center 2-body integral uses a unit function (primitive spherical Gaussian with zero exponent) for function 2. The 2-center 2-body integral also uses the unit function for function 4. These classes of integrals are described by values xx_xx, xs_xx and xs_xs of enum BraKet. To evaluate 2- and 3-center integrals the braket parameter must be provided to the Engine constructor,

BasisSet dfbs("cc-pVDZ-RI", atoms);

Engine eri3_engine(Operator::coulomb, max(obs.max_nprim(), dfbs.max_nprim()),
                   max(obs.max_l(), dfbs.max_l()), 0,
                   std::numeric_limits<real_t>::epsilon(),
                   operator_traits<Operator::coulomb>::default_params(),
                   BraKet::xs_xx);

, or reset using Engine::set method:

Engine eri2_engine(Operator::coulomb, dfbs.max_nprim(), dfbs.max_l());
eri2_engine.set(BraKet::xs_xs);

To compute the integrals provide only the non-unit shells:

// this computes what in chemists notation is denoted (dfbs[0] | obs[1], obs[2])
eri3_engine.compute(dfbs[0], obs[1], obs[2]);

how Engine::compute performs dispatch

Engine::compute exists to simplify the use of the API for beginners, and to allow writing generic code that can handle any operator type, number of shell arguments, and derivative order (e.g. for the purposes of serialization, shell set permutation, etc.). Under the hood it will dispatch the compute command to one of the low-level methods of the Engine object. For the sake of performance, for the 2-body integrals these methods are template functions, with their addresses stored in a table through which Engine::compute performs the dispatch. However, often you can specialize code for a particular operator type, layout, and derivative order; this makes it possible to improve efficiency of dispatch for the 2-body integrals by avoiding the use of Engine::compute and instead calling the underlying compute method for 2-body integrals, Engine::compute2, directly. For example, the eri3_engine.compute(dfbs[0], obs[1], obs[2]) call above can be replaced by eri3_engine.compute2<Operator::coulomb, BraKet::xs_xx, 0>(dfbs[0], Shell::unit(), obs[1], obs[2]). See hartree-fock++.cc for more examples.

shell pair precomputation

As of libint 2.4.0-beta.4 it is possible to provide shell pair data to the Engine object when computing 2-body integrals. This makes possible to avoid the overhead of computing the pair data each time Engine::compute2 is invoked (this typically improves efficiency by few percent). See function compute_2body_fock in hartree-fock++.cc for an example of how to use precomputed shell pair data.

thread-safety

Integral engines are not thread-safe, i.e. their non-const member functions, such as compute() cannot be simultaneously executed from multiple threads. Therefore the user should create one engine per each thread. This topic is outside the scope of this guide. Refer to tests/hartree-fock/hartree-fock++.cc for the examples of how to write parallel code using the libint C++ interface.

Wiki Home

Using Modern C++ API

Contributor Guide

Clone this wiki locally