Integrals#

Exposure to Python layer#

Several integrals are exposed to the Python layer in VeloxChem.

import numpy as np
import veloxchem as vlx

mol_str = """3

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

molecule = vlx.Molecule.from_xyz_string(mol_str)
basis = vlx.MolecularBasis.read(molecule, "cc-pVDZ")

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

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

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()

Electron repulsion integrals#

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

print(g.shape)

(24, 24, 24, 24)

Integral transformation#

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.

# get the MO coefficients
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule, basis)

C = scf_drv.scf_tensors["C_alpha"]

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.

nocc = molecule.number_of_alpha_electrons()

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)