basis_integrals

Compute analytic integrals over contracted Gaussian-type orbitals (GTO).

This module uses pyscf which in turn uses libcint to compute analytic integrals. Basis functions are usually separated into radial and spherical parts,

\[\omega_\mu(r, \theta, \phi) = f_{l_\mu}(r) Y_{l_\mu}^{m_\mu} (\theta, \phi),\]

where \(Y_l^m\) are spherical harmonics. The radial part of contracted Gaussian-type orbitals are given by a monomial times a linear combination of Gaussians,

\[f_l(r) = r^l \sum_{i=1}^N c_i \exp(-\zeta_i r^2).\]

The \(c_i\) are called contraction coefficients and the \(\zeta_i\) scale the width of the Gaussians. \(N\) is the number of primitives. Note that also the nucleus can be modeled using GTOs (nuc_model=1 is point particle).

pyscf molecules have three layers of basis data (see https://pyscf.org/develop/gto_developer.html):

• input, e.g. .atom, .basis

• internal, e.g. ._atom, ._base

• for libcint, e.g. ._atm, ._bas, ._env

The ._env array stores all relevant information about the basis set, including the atom positions, the exponents and contraction coefficients. The properties ._atm and ._bas are used to store the indices of the ._env array. They contain the following information (0 are entries not relevant for us):

• _atm: charge, ind_of_coord, nuc_model, ind_nuc_zeta, 0, 0

• _bas: atom_id, angular_momentum, num_primitives, num_contracted_gtos, 0, ind_exp, ind_contraction_coeff, 0

_append_flat_gaussian(bas, env) → Tuple[ndarray, ndarray]

Append a flat Gaussian to the basis set.

The added flat Gaussian (constant function with value 1) is needed to compute integrals using pyscf.

Parameters:

• bas – The basis array as used by libcint.

• env – The environment array as used by libcint.

Returns:

The modified basis array. env_mod: The modified environment array.

Return type:

bas_mod

_get_one_center_integral(mol: Mole, intor_name: str, components: int = 1) → ndarray

Compute 1-center integrals.

Parameters:

• mol – The molecule containing information about the geometry and the basis set.

• intor_name – The name of the integral to compute. Partial list at pyscf. More options in the c++ library at libcint.

• components – The number of components of the integral. For example, the dipole moment has three components.

Returns:

Vector of integrals.

get_L2_coulomb_matrix(mol: Mole) → ndarray

Compute the two-electron repulsion integral matrix.

Compute

\[D_{\alpha,\beta,\gamma,\delta} = (\eta^\alpha\eta^\beta|\eta^\gamma\eta^\delta)\]

for all basis functions \(\eta^\alpha\), \(\eta^\beta\), \(\eta^\gamma\), and \(\eta^\delta\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Matrix of integrals.(n_basis, n_basis, n_basis, n_basis)

get_basis_dipole(mol: Mole) → ndarray

Compute the basis dipole vector.

Compute

\[d_{i \mu} = \int \omega_\mu(r) r_i \mathrm{d}r\]

for all basis functions \(\omega_\mu\) and coordinates \(r_i\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Matrix of integrals. With shape (3, n_basis).

get_coulomb_matrix(mol: Mole) → ndarray

Compute the 2-center Coulomb integral matrix.

Compute

\[\tilde W_{\mu\nu} = \int\int \frac{\omega_\mu(r_1) \omega_\nu(r_2)}{r_{12}} \mathrm{d} r_1 \mathrm{d} r_2\]

for all basis functions \(\omega_\mu\) and \(\omega_\nu\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Matrix of integrals.

get_coulomb_tensor(mol_rho: Mole, mol_orbital: Mole, shls_slice: tuple | None = None, aosym: str = 's1') → ndarray

Compute the 3-center Coulomb integral tensor between different bases.

Compute the Coulomb integral over one atomic orbital \(\omega_\mu\) of the density basis given by mol_rho and the orbitals \(\omega_\nu\) of the orbital basis given by mol_orbital,

\[\tilde L_{\mu,\alpha\beta} = (\omega_\mu | \eta_\alpha\eta_\beta) = \int\int \frac{\omega_\mu(r_1) \eta_\alpha(r_2) \eta_\beta(r_2)}{r_{12}} \mathrm{d} r_1 \mathrm{d} r_2.\]

Parameters:

• mol_rho – The molecule containing information about the geometry and the density basis set.

• mol_orbital – The molecule containing information about the geometry and the orbital basis set.

Returns:

Tensor of integrals with shape (n_ao_rho, n_ao_orbital, n_ao_orbital).

get_gga_potential_matrix(ao: ndarray, vsigma: ndarray, grid: Grids) → ndarray

Compute the potential matrix for given GGA potential on the grid.

Computes the potential matrix \(V^{(\sigma)}\) from the functional derivative w.r.t. \(\sigma:=\nabla\rho(r)\nabla\rho(r)\). This is necessary for the calculation of the correct gradient of the K/XC functionals. Their gradient can be separated into the usual LDA like gradient and a GGA gradient contribution.

\[\nabla_p E_\mathrm{XC} = \nabla_p E_\mathrm{LDA} + \nabla_p E_\mathrm{GGA}\]

The GGA gradient contribution is given by

\[\nabla_p E_\mathrm{GGA} = \int\underbrace{\frac{\delta E[\rho]}{\delta\sigma(r)}}_{ v_\sigma(r)}\nabla_p \sigma(r)\mathrm{d}r,\]

which can be rewritten as

\[\begin{split}(\nabla_p E_\mathrm{GGA})_\mu &= 2 \sum_\nu p_\nu \underbrace{\int v_\sigma( r)\nabla\omega_\mu(r)\nabla\omega_\nu(r)\mathrm{d}r}_{V^{(\sigma)}_{\mu\nu}}, \\ \nabla_p E_\mathrm{GGA} &= 2 V^{(\sigma)} p.\end{split}\]

This function computes the matrix

\[V^{(\sigma)}_{\mu\nu} = \int v_\sigma(r)\nabla\omega_\mu(r)\nabla\omega_\nu(r)\mathrm{d}r.\]

Parameters:

• ao – The atomic orbitals on the grid with size (n_derivative, n_grid, n_ao).

• vsigma – The functional derivative w.r.t \(\sigma:=\nabla rho \nabla \rho\) on the grid.

• grid – The grid on which the atomic orbitals and the functional are evaluated.

Returns:

Matrix of integrals.

get_normalization_vector(mol: Mole, clip_minimum: float = 1.6e-15) → ndarray

Compute the integral over the basis functions.

Compute

\[w_\mu = \int \omega_\mu(r) \mathrm{d}r\]

for all basis functions \(\omega_\mu\).

Parameters:

• mol – The molecule containing information about the geometry and the basis set.

• clip_minimum – If not None, the integral is clipped from below at this threshold. This should be used to avoid numerical instabilities which can lead to larger errors later on.

Returns:

Vector of integrals.

get_nuclear_attraction_matrix(mol: Mole) → ndarray

Compute the 2-center nuclear attraction integral matrix.

Compute

\[(V_\mathrm{nuc})_{\mu\nu} = \int \eta_\mu(r) \sum_{a=1}^A \frac{Z_a}{\|r - r_a\|} \eta_\nu(r) \mathrm{d} r\]

for all basis functions \(\eta_\mu\) and \(\eta_\nu\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Matrix of integrals.

get_nuclear_attraction_vector(mol: Mole) → ndarray

Compute the nuclear attraction integral vector.

Compute the integral over the potential of the nuclei \(a\)

\[(v_\mathrm{nuc})_\mu = -\int \omega_\mu(r) \sum_{a=1}^A \frac{Z_a}{\|r - r_a\|} \mathrm{d}r\]

for all basis functions \(\omega_\mu\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Vector of integrals.

get_nuclear_gradient_matrix(mol: Mole) → ndarray

Compute the 2-center nuclear gradient integral matrix.

Compute

\[\left(\nabla_{R_A}V_\mathrm{nuc}\right)_{A_{x/y/z}\mu\nu} = \int \omega_\mu(r) \left(\nabla_{R_A}\frac{Z_A}{\|r_1 - A_A\|} \right) \omega_\nu(r) \mathrm{d} r\]

for all basis functions \(\omega_\mu\) and \(\omega_\nu\) and all nuclei \(A\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Matrix of integrals shape (3 * natoms, nao, nao), ordered according to atom id.

get_nuclear_gradient_vector(mol: Mole) → ndarray

Compute the nuclear gradient integral vector.

Compute the integral over the gradient of the potential of the nuclei \(a\)

\[\left(\nabla_{R_A}v_\mathrm{nuc}\right)_{ A_{x/y/z}\mu} = -\int \omega_\mu(r) \nabla_{R_A} \frac{Z_A}{\|r - r_A\|} \mathrm{d}r\]

for all basis functions \(\omega_\mu\) and nuclei \(A\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Vector of integrals, shape (3 * natms, nao), ordered according to atom id.

get_overlap_matrix(mol: Mole) → ndarray

Compute the 2-center overlap matrix.

Compute

\[W_{\mu\nu} = \int \omega_\mu(r) \omega_\nu(r) \mathrm{d}r\]

for all basis functions \(\omega_\mu\) and \(\omega_\nu\).

Parameters:

mol – The molecule containing information about the geometry and the basis set.

Returns:

Matrix of integrals.

get_overlap_tensor(mol_rho: Mole, mol_orbital: Mole) → ndarray

Compute the 3-center overlap integral tensor between different bases.

Compute the overlap integral over one atomic orbital \(\omega_\mu\) of the density basis given by mol_rho and the orbitals \(\eta_\nu\) of the orbital basis given by mol_orbital,

\[\tilde L_{\mu,\alpha\beta} = \langle\omega_\mu | \eta_\alpha\eta_\beta \rangle = \int \omega_\mu(r) \eta_\alpha(r) \eta_\beta(r) \mathrm{d} r.\]

Parameters:

• mol_rho – The molecule containing information about the geometry and the density basis set.

• mol_orbital – The molecule containing information about the geometry and the orbital basis set.

Returns:

Tensor of integrals with shape (n_ao_rho, n_ao_orbital, n_ao_orbital).

get_potential_vector(ao: ndarray, potential: ndarray, grid: Grids) → ndarray

Compute the potential vector for given potential on the grid.

For a given potential \(v_\mathrm{pot}\) and for all atomic orbitals \(\omega_\mu\) compute

\[v_\mu = \sum_{i \in \mathrm{grid}} (\omega_\mu)_i (v_\mathrm{pot})_i w_i \approx \int \omega_\mu(r) v_\mathrm{pot}(r) \mathrm{d}r\]

on the grid. The weights \(w_i\) are taken from grids.

Parameters:

• ao – The atomic orbitals on the grid with size (n_grid, n_ao). If the atomic orbitals also contain derivatives in an array of size (:, n_grid, n_ao), they are ignored.

• potential – The exchange-correlation potential on the grid.

• grid – The grid on which the atomic orbitals and the exchange-correlation potential are defined.

Returns:

Vector of integrals.