Integrals#

Ordering#

In VeloxChem, the atomic orbitals are ordered according to

  1. Orbital angular momentum, \(l\)

  2. Projection of orbital angular momentum, \(m_l\)

  3. The order of atoms in the user input

As an illustration, let us consider a water molecule with an atom ordering and coordinates as follows:

h2o_xyz = """3

O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
"""

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

We employ the cc-pVDZ basis set.

import numpy as np
import veloxchem as vlx

molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "cc-pVDZ", ostream=None)

print(f"Number of basis functions: {basis.get_dimension_of_basis(molecule)}")

Number of basis functions: 24

The basis functions have the following order:

  • \(l=0\), \(m_l=0\):
    \(\chi_{1s}^\mathrm{O}\), \(\chi_{2s}^\mathrm{O}\), \(\chi_{3s}^\mathrm{O}\), \(\chi_{1s}^\mathrm{H_1}\), \(\chi_{2s}^\mathrm{H_1}\), \(\chi_{1s}^\mathrm{H_2}\), \(\chi_{2s}^\mathrm{H_2}\)

  • \(l=1\), \(m_l=-1\):
    \(\chi_{2p_y}^\mathrm{O}\), \(\chi_{3p_y}^\mathrm{O}\), \(\chi_{2p_y}^\mathrm{H_1}\), \(\chi_{2p_y}^\mathrm{H_2}\)

  • \(l=1\), \(m_l=0\):
    \(\chi_{2p_z}^\mathrm{O}\), \(\chi_{3p_z}^\mathrm{O}\), \(\chi_{2p_z}^\mathrm{H_1}\), \(\chi_{2p_z}^\mathrm{H_2}\)

  • \(l=1\), \(m_l=+1\):
    \(\chi_{2p_x}^\mathrm{O}\), \(\chi_{3p_x}^\mathrm{O}\), \(\chi_{2p_x}^\mathrm{H_1}\), \(\chi_{2p_x}^\mathrm{H_2}\)

  • \(l=2\), \(m_l=\{-2, -1, 0, +1, +2\}\):
    \(\chi_{3d_{xy}}^\mathrm{O}\), \(\chi_{3d_{yz}}^\mathrm{O}\), \(\chi_{3d_{z^2}}^\mathrm{O}\), \(\chi_{3d_{xz}}^\mathrm{O}\), \(\chi_{3d_{x^2 - y^2}}^\mathrm{O}\)

Exposure to Python layer#

Several integrals are exposed to the Python layer in VeloxChem.

Overlap#

overlap_drv = vlx.OverlapIntegralsDriver()
S = overlap_drv.compute(molecule, basis).to_numpy()

Kinetic energy#

kinetic_drv = vlx.KineticEnergyIntegralsDriver()
T = kinetic_drv.compute(molecule, basis).to_numpy()

Nuclear potential#

nucpot_drv = vlx.NuclearPotentialIntegralsDriver()
V = -1.0 * nucpot_drv.compute(molecule, basis).to_numpy()

One-electron Hamiltonian#

h = T + V

Linear momentum#

linmom_drv = vlx.LinearMomentumIntegralsDriver()
linmom_mats = linmom_drv.compute(molecule, basis)

p_x = linmom_mats.x_to_numpy()
p_y = linmom_mats.y_to_numpy()
p_z = linmom_mats.z_to_numpy()

Electric dipole moment#

dipole_drv = vlx.ElectricDipoleIntegralsDriver()
dipole_mats = dipole_drv.compute(molecule, basis)

mu_x = -1.0 * dipole_mats.x_to_numpy()
mu_y = -1.0 * dipole_mats.y_to_numpy()
mu_z = -1.0 * dipole_mats.z_to_numpy()

Magnetic dipole moment#

angmom_drv = vlx.AngularMomentumIntegralsDriver()
angmom_mats = angmom_drv.compute(molecule, basis)

m_x = -0.5 * angmom_mats.x_to_numpy()
m_y = -0.5 * angmom_mats.y_to_numpy()
m_z = -0.5 * angmom_mats.z_to_numpy()

Electron repulsion integrals#

eri_drv = vlx.ElectronRepulsionIntegralsDriver()
g = eri_drv.compute_in_memory(molecule, basis)

print(g.shape)

(24, 24, 24, 24)

Integral transformations#

# number of occupied orbitals
nocc = molecule.number_of_alpha_electrons()

# get the MO coefficients
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()

scf_results = scf_drv.compute(molecule, basis)

C = scf_results["C_alpha"]

One-electron integrals#

The transformation of one-electron integrals from the AO to the MO basis takes the form

\[ \Omega_{pq}^\mathrm{MO} = \sum_{\alpha, \beta} c^*_{\alpha p} \Omega_{\alpha\beta}^\mathrm{AO} c_{\beta q} \]
mu_z_mo = np.einsum("ap, ab, bq -> pq", C, mu_z, C)

print(f"Electronic dipole moment: {2 * mu_z_mo[:nocc,:nocc].trace():10.6f} a.u.")

Electronic dipole moment:  -1.391848 a.u.

Two-electron integrals#

Most correlated wave function methods rely on integrals in the molecular orbital basis. For some methods, like Møller–Plesset second order perturbation theory, this step is the most time-consuming step of the calculation.

Conventional#

The transformation of the AO integrals to the MO basis can be done as follows

\[\begin{eqnarray*} ( \alpha \beta | \gamma s ) &= \sum_{\delta} c_{\delta s} ( \alpha \beta | \gamma \delta ) \\ % ( \alpha \beta | r s ) &= \sum_{\gamma} c^*_{\gamma r} ( \alpha \beta | \gamma s ) \\ % ( \alpha q | r s ) &= \sum_{\beta} c_{\beta q} ( \alpha \beta | r s ) \\ % ( p q | r s ) &= \sum_{\alpha} c^*_{\alpha p} ( \alpha q | r s ) \\ \end{eqnarray*}\]

The computational cost of this procedure is \(O(N^5)\), since each summation involves five indices. Note that the intermediate result of the transformation needs to be explicitly stored in memory. This can be demanding as the required memory increases as \(O(N^4)\), where \(N\) is the number of contracted basis functions.

g_1 = np.einsum("ds, abcd -> abcs", C, g)
g_2 = np.einsum("cr, abcs -> abrs", C, g_1)
g_3 = np.einsum("bq, abrs -> aqrs", C, g_2)

g_mo = np.einsum("ap, aqrs -> pqrs", C, g_3)

print(g_mo.shape)

(24, 24, 24, 24)

Specific blocks a ERIs are made available directly in the MO basis, both in the physicist’s and the chemist’s notation.

erimo_drv = vlx.MOIntegralsDriver()

phys_oovv = erimo_drv.compute_in_memory(molecule, basis, scf_drv.mol_orbs, "phys_oovv")
print("[oo|vv]:", phys_oovv.shape)

chem_ovov = erimo_drv.compute_in_memory(molecule, basis, scf_drv.mol_orbs, "chem_ovov")
print("(ov|ov):", chem_ovov.shape)

np.testing.assert_allclose(phys_oovv, chem_ovov.swapaxes(1, 2), atol=1e-15)

[oo|vv]: (5, 5, 19, 19)

(ov|ov): (5, 19, 5, 19)

Fock-matrix driven#

An alternative way of getting the molecular orbital integrals is through the direct formation of many Fock matrices. Taking the ovov-block of MO integrals as an example, we can in practice build \(N_{occ} \times N_{occ}\) density matrices using the coefficients of the occupied orbitals, and form the corresponding Fock matrices that can then be transformed into molecular basis by the coefficients of the virtual orbitals.

\[\begin{eqnarray*} D^{ij}_{\alpha \beta} &=& c^*_{\alpha i} c^*_{\beta j} \\ % K^{ij}_{\gamma \delta} &=& \sum_{\alpha, \beta} D^{ij}_{\alpha\beta} (\alpha \gamma| \beta \delta ) \\ % (i a|j b ) &=& \sum_{\gamma, \delta} c_{\gamma a} K^{ij}_{\gamma \delta} c_{\delta b} \end{eqnarray*}\]

The computational cost of this approach is formally \(O(N^6)\); however, in practice the cost scales between \(O(N^4)\) and \(O(N^5)\) due to screening of integrals in the formation of Fock matrices. An advantage of this approach is that the Fock matrices can be computed and stored on individual compute nodes, making it suitable for large-scale parallelization on HPC systems.

VeloxChem uses this approach to serve electron repulsion integrals (ERIs) for the correlated wave function programs Gator and MultiPsi.

i, j = 0, 3
D_ij = np.einsum("a, b -> ab", C[:, i], C[:, j])
K_ij = np.einsum("acbd, ab -> cd", g, D_ij)
g_ivjv = np.einsum("ca, cd, db -> ab", C[:, nocc:], K_ij, C[:, nocc:])

print(g_ivjv.shape)
np.testing.assert_allclose(g_ivjv, chem_ovov[i, :, j, :], atol=1e-15)

(19, 19)