Integrals#

import numpy as np
import veloxchem as vlx

Ordering#

In VeloxChem, the atomic orbitals are ordered according to

  1. Orbital angular momentum, l

  2. Projection of orbital angular momentum, ml

  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
"""

We employ the cc-pVDZ basis set.

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()}")
Number of basis functions: 24

The basis functions have the following order:

  • l=0, ml=0:
    χ1sO, χ2sO, χ3sO, χ1sH1, χ2sH1, χ1sH2, χ2sH2

  • l=1, ml=1:
    χ2pyO, χ3pyO, χ2pyH1, χ2pyH2

  • l=1, ml=0:
    χ2pzO, χ3pzO, χ2pzH1, χ2pzH2

  • l=1, ml=+1:
    χ2pxO, χ3pxO, χ2pxH1, χ2pxH2

  • l=2, ml={2,1,0,+1,+2}:
    χ3dxyO, χ3dyzO, χ3dz2O, χ3dxzO, χ3dx2y2O

Exposure to Python layer#

Several integrals are exposed to the Python layer in VeloxChem.

Overlap#

S = vlx.compute_overlap_integrals(molecule, basis)

Kinetic energy#

T = vlx.compute_kinetic_energy_integrals(molecule, basis)

Nuclear potential#

V = vlx.compute_nuclear_potential_integrals(molecule, basis)

One-electron Hamiltonian#

h = T + V

Linear momentum#

linmom_mats = vlx.compute_linear_momentum_integrals(molecule, basis)

p_x = linmom_mats[0]
p_y = linmom_mats[1]
p_z = linmom_mats[2]

Electric dipole moment#

dipole_mats = vlx.compute_electric_dipole_integrals(molecule, basis)

mu_x = dipole_mats[0]
mu_y = dipole_mats[1]
mu_z = dipole_mats[2]

Magnetic dipole moment#

angmom_mats = vlx.compute_angular_momentum_integrals(molecule, basis)

m_x = -0.5 * angmom_mats[0]
m_y = -0.5 * angmom_mats[1]
m_z = -0.5 * angmom_mats[2]

Electron repulsion integrals#

fock_drv = vlx.FockDriver()
g = fock_drv.compute_eri(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#

AO to MO transformation#

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

ΩpqMO=α,βcαpΩαβAOcβq

or, equivalently,

ΩMO=CΩAOC
mu_z_mo = np.einsum("ap, ab, bq -> pq", C, mu_z, C)

MO to AO transformation#

To do the reversed operation, we make use of the relation

CSC=I

from which we deduce

C1=CS[C]1=SC

and we get

ΩAO=SCΩMOCS
mu_z_ao = np.einsum("ab, bp, pq, cq, cd -> ad", S, C, mu_z_mo, C, S)

Test the forth-and-back transformations against the original AO matrix.

np.testing.assert_allclose(mu_z_ao, mu_z, atol=1e-12)

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, the step of transforming integrals from the AO to the MO basis can the most time-consuming step of the calculation.

Conventional#

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

(αβ|γs)=δcδs(αβ|γδ)(αβ|rs)=γcγr(αβ|γs)(αq|rs)=βcβq(αβ|rs)(pq|rs)=αcαp(αq|rs)

The computational cost of this procedure is O(N5), 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(N4), 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 Nocc×Nocc 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.

Dαβij=cαicβjKγδij=α,βDαβij(αγ|βδ)(ia|jb)=γ,δcγaKγδijcδb

The computational cost of this approach is formally O(N6); however, in practice the cost scales between O(N4) and O(N5) 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)