using modern CPlusPlus API (legacy)

=============================================================================

intro

libint library provides 2 types of interfaces:

1. low-level, i.e. intended for use from any native language, including Fortran77. This is intended for experts who want to integrate libint library's advanced features (like vectorization) and/or want to squeeze the last drop of performance from libint.
2. high-level based on "modern C++" (C++11 or higher); this API is intended for use by everyone except the few obsessed souls. It is still very efficient, but some optimizations are simply not possible due to the tension between performance and easy-to-use interface. It also does not expose ALL functionality of libint at the moment.

The easiest (and best for 99.9% of users) way to get started with libint is to use its "modern C++" interface. Even if you are not a C++ programmer, this is the recommended way forward.

This document describes the features of the C++ libint interface. 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 parallel closed-shell Hartree-Fock energy (with efficiency comparable to production codes): tests/hartree-fock/hartree-fock++.cc

=============================================================================

prerequisites

1. C++ compiler that supports C++11 standard. C++11 standard is the second most recent C++ international standard for C++, hence most modern compilers support it. You will need to use -std=c++11 or -std=c++0x compiler flag to turn this support on.
2. libint library, configured (./configure ...), compiled (make) and installed (make install). To get started you could try this library . You will need to specify 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 init(). 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::init(); // safe to use libint now

  // all other code snippets go here

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

  return 0;
}

P.S. It's a good idea to wrap these calls in a singleton class.

create basis set

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

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

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

Think of libint2::BasisSet as a slightly decorated std::vector<libint2::Shell>. This means that you can manipulate it as a std::vector directly, e.g. append shells, reorder shells, sort shells, etc. libint2::BasisSet exists to make basis set construction easier (e.g. by reading standard basis sets from a stash accompanying libint). However it is very easy to construct shells and basis sets directly, as shown in tests/hartree-fock/hartree-fock.cc .

libint2::BasisSet uses 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:

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

or even simpler:

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

create integral engine

Here's how to create a one-body integral engine for computing overlap integrals:

// obs was created above

libint2::OneBodyEngine onebody_engine(libint2::OneBodyEngine::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 OneBodyEngine is a value of OneBodyEngine::operator_type type that specifies the operator (or operator set; see below) that the engine will evaluate. The currently supported values are OneBodyEngine::overlap (overlap integrals), OneBodyEngine::kinetic (kinetic energy), OneBodyEngine::nuclear (Coulomb potential due to point charges), OneBodyEngine::emultipole1 (overlap + electric dipole), OneBodyEngine::emultipole2 (operators in emultipole1 + Cartesian electric quadrupole), and OneBodyEngine::emultipole3 (operators in emultipole2 + Cartesian electric octupole).

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.

The fourth parameter to the constructor (not shown) is the order of geometrical derivatives. This feature is not ready for production use yet.

Making engines for some types of one-body 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:

libint2::OneBodyEngine onebody_engine2(libint2::OneBodyEngine::nuclear, // will compute nuclear attraction ints
                                       obs.max_nprim(), obs.max_l());
onebody_engine.set_params(atoms); // atoms specifies the charge and position of each nucleus

Another example are electric (Cartesian) multipole integrals. The default is to use the 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 std::array<double,3>, via set_params:

libint2::OneBodyEngine emultipole1_engine(libint2::OneBodyEngine::emultipole1, // will compute overlap + electric dipole moments
                                          obs.max_nprim(), obs.max_l());
emultipole1_engine.set_params(std::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, except the engine type is a template class parametrized by the integral kernel type (this may change before the API is finalized). For example, two-electron repulsion ints engine:

typedef libint2::TwoBodyEngine<libint2::Coulomb> coulomb_engine_t;
coulomb_engine_t twobody_engine(obs.max_nprim(), obs.max_l() );

Actually, the constructor for TwoBodyEngine class takes additional parameters, but all of them have default values:

typedef libint2::TwoBodyEngine<libint2::Coulomb> coulomb_engine_t;
coulomb_engine_t twobody_engine(obs.max_nprim(), obs.max_l(),
                                0,  // order of geometrical derivatives;
                                    // 0 = compute regular integrals (default)
                                    // 1 = compute first-order derivative integrals
                                    // only regular derivatives are supported at the moment
                                std::numeric_limits<double>::epsilon() // target precision of the integrals
                                                                       // default = as close to machine precision as possible
                               );

Several other kernel types are supported, all of these appear in explicitly correlated methods. For example, to compute integrals over a linear combination of Gaussian-Type Geminals (GTG), f(r12) = 0.1 exp(-0.2 r12) - 0.3 exp(-0.4 r12) + 0.5 exp(-0.6 r12), use this code:

typedef libint2::TwoBodyEngine<libint2::cGTG> contracted_gtg_engine_t;
std::vector<std::pair<double,double>> cgtg_params{{0.2,0.1}, {0.4,-0.3}, {0.6, 0.5}};
contracted_gtg_engine_t twobody_engine2(obs.max_nprim(), obs.max_l(),
                                        0, std::numeric_limits<double>::epsilon(),
                                        cgtg_params // geminal parameters
                                       );

Other kernel types are used similarly:

typedef libint2::TwoBodyEngine<libint2::cGTG_times_Coulomb> cgtg_coulomb_engine_t;
cgtg_coulomb_engine_t twobody_engine3(obs.max_nprim(), obs.max_l(),
                                      0, std::numeric_limits<double>::epsilon(),
                                      cgtg_params);
typedef libint2::TwoBodyEngine<libint2::DelcGTG_square> delcgtg2_engine_t;
delcgtg2_engine_t twobody_engine4(obs.max_nprim(), obs.max_l(),
                                  0, std::numeric_limits<double>::epsilon(),
                                  cgtg_params);

compute integrals

Using the engines to compute integrals is trivial. Here's how to compute the entire set of one electron 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
                                // ...

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

    std::cout << "compute shell set {" << s1 << "," << s2 << "} ... ";
    const auto* ints_shellset = onebody_engine2.compute(obs[s1], obs[s2]);
    std::cout << "done" << std::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

    // 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)
        std::cout << "  " << bf1+f1 << " " << bf2+f2 << " " << ints_shellset[f1*n2+f2] << std::endl;
  }
}

Each call to OneBodyEngine::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 emultipole1_engine will return multiple (namely, four) consecutive shell blocks, one block per corresponding operator in the integral set computed by that engine. The 4 shell blocks correspond to (in order) overlap, x-dipole, y-dipole, and z-dipole operators. Here's a code snippet that shows how to digest the integrals:

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

    std::cout << "compute shell set {" << s1 << "," << s2 << "} ... ";
    const auto* ints_shellset = emultipole1_engine.compute(obs[s1], obs[s2]);
    std::cout << "done" << std::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

    const auto* overlap_shellset = ints_shellset + n1*n2;
    const auto* edipole_x_shellset = overlap_shellset + n1*n2;
    const auto* edipole_y_shellset = edipole_x_shellset + n1*n2;
    const auto* edipole_z_shellset = edipole_y_shellset + n1*n2;
    std::cout << "e.g. z dipole integrals:" << std::endl;
    for(auto f1=0; f1!=n1; ++f1)
      for(auto f2=0; f2!=n2; ++f2)
        std::cout << "  " << bf1+f1 << " " << bf2+f2 << " " << edipole_z_shellset[f1*n2+f2] << std::endl;
  }
}

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

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

TwoBodyEngine objects are used similarly. More notes to come soon.

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.