Integrals#
import numpy as np
import veloxchem as vlx
Ordering#
In VeloxChem, the atomic orbitals are ordered according to
Orbital angular momentum,
Projection of orbital angular momentum,
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:
, :
, , , , , , , :
, , , , :
, , , , :
, , , , :
, , , ,
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
or, equivalently,
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
from which we deduce
and we get
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
The computational cost of this procedure is
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
The computational cost of this approach is formally
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)