There are two main ways to use the libecpint library. In the high-level API, you pass details of the system (basis set, coordinates, ECPs) to libecpint, and it handles the computation of all the integrals and/or derivatives automatically, returning arrays of the (Cartesian) integrals. In the low-level API, you control the calculation of the integrals yourself, calling the relevant routines as and when you need them. We envisage that for most purposes, the high-level API will be more appropriate, and is easier to use.

Important note about ECP definitions

ECP powers (u_ns) are assumed to be including the r^2 from the spherical Jacobian, as is the convention followed with the Stuttgart-Dresden ECPs. In the code, we subtract two from this (see the ECP constructor) as integration is done in Cartesian coordinates. In general, the inputted n should only ever be 0, 1, or (most commonly) 2. If your ECP definitions are instead -2 or -1, or predominantly 0, that suggests you are following the other convention of not including the Jacobian, so you should add 2. Yes, this is annoying; if in doubt check the built-in ECP library for examples.

High Level API

Examples of using this API can be found in tests/lib/api_test1 and tests/lib/api_test2.

The ECPIntegrator Object

The first step in using this API is to include the libecpint header and create an ECPIntegrator object:

include <libecpint.hpp>
ECPIntegrator factory;

This object will form the main interface to all of the subsequent routines, and is described in detail in the Library API.


There are three steps to initialising the ECPIntegrator before it can be used to calculate integrals. These are:

  1. specifying the Gaussian basis set

  2. specifying the ECP basis

  3. calling the init routine

These steps are performed as follows:

factory.set_gaussian_basis(N_shells, g_coords, g_exps, g_coefs, g_ams, g_lengths);
factory.set_ecp_basis(N_ecps, u_coords, u_exps, u_coefs, u_ams, u_ns, u_lengths);

where N_shells and N_ecps are the numbers of shells in the Gaussian basis, and the number of ECP centers, respectively, while deriv_order is the maximum derivatives needed (0, 1, or 2). The rest of the parameters are:

  • the Cartesian coordinates in Bohr (g_coords, u_coords);

  • the exponents, coefficients, and angular momenta (and powers, u_ns, for the ECPs) comprising the basis sets

  • the number of exponents per shell (g_lengths, u_lengths)

tests/lib/api_test1 shows how to specify these for HBr in the AVDZ/AVDZ-PP basis.

NOTE: The atom order given in set_gaussian_basis fixes the atom order for all the derivatives, as will be described later.

The ECP library

Alternatively, step 2 can be replaced by reading the ECPs from the built-in library provided with libecpint. This can be found in share/libecpint. To do this, you call:

factory.set_ecp_basis_from_library(N_ecps, u_coords, u_charges, u_names, share_dir);

The new parameters are:

  • u_charges a list of atomic numbers for the ECPs, corresponding to the centers in u_coords;

  • u_names the ECP names for each ECP, e.g. ecp10mdf;

  • share_dir the absolute path to the share/libecpint directory, which must be passed by you.

The currently available ECPs in the library (more being added soon), and the atoms they are available for, are given below:




K – Kr (Z = 19 – 36)


Rb – Xe (Z = 37 – 54)


Cs, Ba (Z = 55, 56)


Hf – Rn (Z = 72 – 86), Ac – U (Z = 89 – 92)


Fr, Ra (Z = 87, 88)


Na – La (Z = 11 – 57), Hf – Bi (Z = 72 – 83), U – Pu (Z = 92 – 94)


  1. Put the ECP in MOLPRO format in share/libecpint/raw as NAME.ecp, where NAME is the name of the ECP; make sure that any exponents are with E (C-convention) not D (Fortran-convention).

  2. Make the top line of NAME.ecp be the NAME.

  3. In share/libecpint run python3 NAME (Python >=3.6 required, with lxml module). This will create NAME.xml in the share/libecpint/xml folder, and this ECP will now be available for use by libecpint. Please consider creating a pull request so that everyone can benefit from the addition!

Note that the n value for each ECP primitive should typically be 0, 1, or 2 (for the Stuttgart-Dresden ECPs, for example, it is always 2). Some input formats follow a convention of subtracting or adding 2 to this.

Computing integrals

Computing integrals over all shell pairs is then very simple:


To retrieve these you then create a shared pointer to a vector:

std::shared_ptr<std::vector<double>> integrals = factory.get_integrals();
double I00 = (*ints)[0]; // example for accessing element (0, 0)

These are stored in row-order, and are in a Cartesian Gaussian basis. Typically these would be converted to a spherical harmonic Gaussian basis (we might add the ability to do this later). We follow canonical Cartesian order, so for a d-type function this would be xx, xy, xz, yy, yz, zz, and the order of the shells is the same as when you called set_gaussian_basis. The total number of Cartesian gaussians is stored in factory.ncart; you can access the ij-th integral as

double Iij = (*ints)[i*factory.ncart + j]

First derivatives

First derivatives are similarly calculated by calling


Note that this will only work if init was called with deriv_order > 0. This will return an array of 3*factory.natoms shared pointers to the integral derivatives with respect to each coordinate. The order is x, y, z, and the order of atoms matches that specified in set_gaussian_basis. For example, to get the array of integral derivatives with respect to the y-coordinate of the n-th atom, you would do:

std::vector<std::shared_ptr<std::vector<double>>> first_derivs = factory.get_first_derivs();
I_dy_atom_n_00 = (*first_derivs[3*n+1])[0];

The order of the elements in each array is identical to that from compute_integrals.

Second derivatives

As for first derivatives, second derivatives are computed as


and are provided as a vector of shared pointers to arrays. The order of these derivatives is somewhat more complicated though, and takes full advantage of symmetry. If the atoms are A, B, C, … as specified in set_gaussian_basis, then they are blocked as follows:

AA, AB, AC, ..., BB, BC, ..., CC, ...

Within each block the order is xx, xy, xz, yy, yz, z on the diagonal (e.g. AA, BB, CC, …) and xx, xy, xz, yx, yy, yz, zx, zy, zz on the off-diagonal (e.g. AB, AC, BC, …). There is a helper macro for this, H_MACRO, defined in api.hpp. So for example to get the derivative with respect to Ay (atom index 0) and Cx (atom index = 2) in a system with N atoms, you would do

int deriv_index = H_START(0, 2, N) + 3; // yx is index 3 for off-diagonal blocks
std::shared_ptr<std::vector<double>> h_Ay_Cx = factory.get_second_derivs()[deriv_index];

Hopefully this doesn’t give you too much of a headache working out.

Updating coordinates

To update the coordinates for the basis and ECPs (for example, after a step in a geometry optimisation), simply pass the new coordinates in the same order they were given when initialised. This is done as:

factory.update_gaussian_basis_coords(N_shells, g_coords);
factory.update_ecp_basis_coords(N_ecps, u_coords);

You then call the compute routines when you need the new integrals and/or derivatives.

NOTE you will need to re-get the pointers using the get routines every time you recompute integrals/derivatives.


TODO detail optional settings that can be passed to ECPIntegrator

Low Level API

Examples of using this API can be found in tests/lib/[name]_test[number] where [name] is int (integrals), deriv (first derivatives), or hess (second derivatives), and [number] is 1 or 2.

The ECPIntegral Object

To be able to calculate integrals and derivatives at the shell-pair level, you need to create an ECPIntegral object instead. This is done as follows:

#include <libecpint/ecpint.hpp>
ECPIntegral ecpint(maxLB, maxLU, deriv_order);

where maxLB is the maximum angular momentum in the Gaussian basis, maxLU is the maximum angular momentum in the ECP basis, and deriv_order is the maximum order of derivative needed (defaults to 0).

Making shells and ECPs

The compute functions in this API require ECP and GaussianShell objects representing the ECP and Gaussian basis functions, respectively. These are populated for the ECP as follows:

#include <libecpint/ecp.hpp>
double C[3] = {Cx, Cy, Cz};
ECP newU(C);
newU.addPrimitive(n, l, x, c);
// addPrimitive for each primitive in ECP

where C is the coordinates of the center of the ECP (in Bohr), and n, l, x, c are the power, angular momentum, exponent, and coefficient of each primitive in that ECP. For the Gaussian basis instead:

#include <libecpint/gshell.hpp>
double A[3] = {Ax, Ay, Az};
GaussianShell shellA(A, l);
shellA.addPrim(x, c);
// addPrim for each primitive in this shell

where A is the coordinates of the center of the shell (in Bohr), and l, x, c are as above. You need to either create these objects for every shell and ECP as they are needed, or store them, so that they can be passed to the compute routines below.

Computing integrals

Computing integrals over a shell pair is then fairly simple. You call

#include <libecpint/mathutil.hpp>
TwoIndex<double> results;
ecpint.compute_shell_pair(newU, shellA, shellB, results);

This will place the integrals in a matrix, results, with dimensions (ncartA, ncartB). The order in each is the canonical Cartesian order as described earlier for the high-level API. Elements of the results matrix can be accessed in two ways:

double Iij = results(i, j);
// or
double Iij =[i*ncartB + j];

First derivatives

Similarly, first derivatives are calculated using

std::vector<TwoIndex<double>> results;
ecpint.compute_shell_pair_derivative(newU, shellA, shellB, results);

Now results is an array of 9 matrices (the matrices are ordered the same as for compute_shell_pair above). These are derivative matrices with respect to Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz, where A, B, C are the centers of shellA, shellB, and newU respectively.

NOTE These derivatives are designed to be additive. Thus, if A==B, then the derivative for the Ax coordinate will be the sum Ax+Bx, etc.

Second derivatives

The second derivatives are calculated in the same simple manner, but their ordering is, as with the high-level API, much more complicated.

std::vector<TwoIndex<double>> results;
ecpint.compute_shell_pair_second_derivative(newU, shellA, shellB, results);

Now results contains 45 different derivative matrices. These are, in order:

AxAx, AxAy, AxAz, AyAy, AyAz, AzAz,
AxBx, AxBy, AxBz, AyBx, AyBy, AyBz, AzBx, AzBy, AzBz,
AxCx, AxCy, AxCz, AyCx, AyCy, AyCz, AzCx, AzCy, AzCz,
BxBx, BxBy, BxBz, ByBy, ByBz, BzBz,
BxCx, BxCy, BxCz, ByCx, ByCy, ByCz, BzCx, BzCy, BzCz,
CxCx, CxCy, CxCz, CyCy, CyCz, CzCz

where A, B, C are as described earlier. These are again additive but THERE IS NOW SYMMETRY TO CONSIDER. I strongly recommend that you look at the code in api.cpp for the construction of the full Hessian, to see how this symmetry has to be taken account of when any of the three centers are the same, as it is too complicated to describe here.


TODO detail optional settings that can be passed to ECPIntegral