Molecular Integrals
1. Hartree-Fock Theory
2. Second Quantization
3. Many-Body Perturbation Theory
4. Configuration Interaction with Single Excitations
5. Static Response Properties
6. Molecular Integrals
Molecular Integrals¶
The Jupyter notebook for this tutorial can be found here.
Electronic structure methods for molecular calculations are usually formulated within a basis of atom-centered Gaussian-type orbitals (GTOs). Integrals involving one- and two-body quantum mechanical operators and GTOs are an essential component of these methods. This project demonstrates the use of recursion relationships to evaluate integrals over simple one-body operators. Python code for evaluating overlap, dipole, quadrupole, linear momentum, angular momentum, and kinetic energy integrals are provided. This tutorial closely follows the first few sections of Chapter 9 of Molecular Electronic-Structure Theory by Helgaker, Jørgensen, and Olsen.
Contracted Gaussian-Type Orbitals¶
We consider contracted GTOs that are formed by linear combinations of primitive Gaussian orbitals centered at different atomic positions, ${\bf A}$: $$\begin{align} G_{ijk}({\bf r}, a, {\bf A}) = x^i_A y_A^j z_A^k \exp(-a r_A^2) \end{align}$$ Here, ${\bf r}$ represents the electronic coordinate, and ${\bf r}_A = {\bf r} - {\bf A}$. The superscripts $i$, $j$, and $k$ define the angular momentum of the function, $l = i + j + k$, and the positive, real-valued exponent factor $a$ defines the tightness of the GTO. In this tutorial, we limit our considerations to Cartesian GTOs, where there are $N = \frac{(l+1)(l+2)}{2}$ GTOs in a shell with total angular momentum $l$ (1 GTO for $s$-type functions, 3 for $p$-type functions, 6 for $d$-type functions, etc.). Integrals over Cartesian GTOs can be transformed into a spherical GTO basis after they have been evaluated.
A nice property of the Cartesian GTOs is that they can be factorized as $$\begin{align} G_{ijk}({\bf r}, a, {\bf A}) = G_{i}(x, a, A_x) G_{j}(y, a, A_y) G_{k}(z, a, A_z) \end{align}$$ with $$\begin{align} G_i(x, a, A_x) = x_A^i \exp(-ax_A^2) \end{align}$$ etc. As a result, integrals over Cartesian GTOs are separable and can be carried out considering one Cartesian component at a time.
Obtaining Contraction Coefficients and Exponents¶
Files containing the contraction coefficients and exponents defining the GTOs for standard basis sets can be found at the Basis Set Exchange. If you visit the Basis Set Exchange, you will find that the files containing the basis set information can be formatted in multiple ways. For example, if we select the Dunning-type correlation consistent polarized valence double-$\zeta$ basis set (cc-pVDZ) for a carbon atom, in Psi4 format, the basis set file looks like:
spherical
!----------------------------------------------------------------------
! Basis Set Exchange
! Version 0.12
! https://www.basissetexchange.org
!----------------------------------------------------------------------
! Basis set: cc-pVDZ
! Description: cc-pVDZ
! Role: orbital
! Version: 1 (Data from ccRepo/Grant Hill)
!----------------------------------------------------------------------
****
C 0
S 9 1.00
6.665000D+03 6.920000D-04
1.000000D+03 5.329000D-03
2.280000D+02 2.707700D-02
6.471000D+01 1.017180D-01
2.106000D+01 2.747400D-01
7.495000D+00 4.485640D-01
2.797000D+00 2.850740D-01
5.215000D-01 1.520400D-02
1.596000D-01 -3.191000D-03
S 9 1.00
6.665000D+03 -1.460000D-04
1.000000D+03 -1.154000D-03
2.280000D+02 -5.725000D-03
6.471000D+01 -2.331200D-02
2.106000D+01 -6.395500D-02
7.495000D+00 -1.499810D-01
2.797000D+00 -1.272620D-01
5.215000D-01 5.445290D-01
1.596000D-01 5.804960D-01
S 1 1.00
1.596000D-01 1.000000D+00
P 4 1.00
9.439000D+00 3.810900D-02
2.002000D+00 2.094800D-01
5.456000D-01 5.085570D-01
1.517000D-01 4.688420D-01
P 1 1.00
1.517000D-01 1.000000D+00
D 1 1.00
5.500000D-01 1.0000000
****
At the top line, "spherical" indicates that the basis was optimized for spherical GTOs, as opposed to Cartesian GTOs. Further down, we see the atom name, "C", followed by sets of data labeled with different symbols for the angular momentum quantum number. The cc-pVDZ basis is a "split-valence" type basis, where there is one contracted GTO for the core (1$s$) orbital and two contracted GTOs for each of the valence orbital subshells (2$s$ and 2$p$). The first "S" label corresponds to the core function; we can see that this function is defined by a linear combination of 9 primitive GTO functions, with the exponents and contraction coefficients defined by the first and second columns. Next, there are coefficents and exponents defining the primitive GTOs for two $s$-type functions and two $p$-type functions that represent the valence orbitals. Note that, even though a neutral carbon atom only has 4 $p$-type electrons, the full set of $p$-type functions is considered as the valence space and share the same contraction coefficients and exponents. Lastly, we see a single additional $d$-type function, which is called a "polarization" function. Polarization functions are characterized by higher angular momentum higher than the valence functions.
In principle, we could download basis set files for all of the atoms we want to consider and write code to parse the basis set files and store the contraction coefficients and exponents for each atom. However, that activity is not particularly exciting and is not the focus of this tutorial. Instead, we will use the existing machinery in the Psi4 package to parse the basis set for us. The following code snippet imports the Psi4 Python library, defines the geometry for a water molecule, defines a basis set (in this case, cc-pVDZ), and creates the corresponding "wavefunction" and "basis" objects. The basis object will contain all of the information we need to evaluate our molecular integrals.
import psi4
mol = psi4.geometry("""
0 1
O
H 1 1.0
H 1 1.0 2 104.5
symmetry c1
""")
# activate the molecule
psi4.activate(mol)
# set the options
psi4.set_options({'basis': 'cc-pvdz',
'puream': 'False',
})
# suppress psi4 printing
psi4.core.be_quiet()
# we can use the wavefunction object to process the basis set for us
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('BASIS'))
# basis set object
basis = wfn.basisset()
# the number of basis functions
nbf = basis.nbf()
# the number of shells
nshell = basis.nshell()
Note that we have added the option 'puream': 'False' to the Psi4 options dictionary. This option ensures that Psi4 uses Cartesian GTOs, even though the cc-pVDZ basis is optimized for spherical GTOs. In this way, we will be able to compare our integrals over Cartesian GTOs to the integrals evaluated by the Psi4 package.
Gaussian Product Theorem¶
The product of two Gaussian functions at centers A and B is given by another Gaussian function at a different center, P. As an example, let us consider the product of the $x$-component of two GTOs of the form defined above $$\begin{align} G_i(x, a, A_x)G_j(x, b, B_x) &= x_A^i x_B^j \exp(-ax_A^2) \exp(-b x_B^2) \\ &= x_A^i x_B^j K_{ab}^x \exp(-px_P^2) \end{align}$$ where $p = a + b$ is the total exponent, $P_x$ represents a center-of-charge coordinate $$\begin{align} P_x = \frac{a A_x + b B_x}{p} \end{align}$$ $X_{AB}$ is the relative coordinate defined by $$\begin{align} X_{AB} = A_x - P_x \end{align}$$ and $K_{ab}^x$ is the pre-exponential factor $$\begin{align} K_{ab}^x = \exp(-\mu X_{AB}^2) \end{align}$$ The symbol $\mu = \frac{ab}{a+b}$ in the pre-exponential factor is the reduced exponent.
One-Electron Integral Evaluation with the Obara-Saika Scheme¶
The Obara-Saika scheme defines molecular integrals with recursion relations derived from the translational invariance of the integrals. A detailed derivation of the recursion relations for various one-electron integrals can be found in section 9.3 of Helgaker's book.
Overlap Integrals¶
We start with the simplest molecular integral, which is the overlap of two primitive Cartesian Gaussian functions, $\langle G_{a} | G_{b} \rangle$, where we have introduced short-hand notation $$\begin{align} G_a({\bf r}) &= G_{ikm}({\bf r}, a, {\bf A})\\ G_b({\bf r}) &= G_{jln}({\bf r}, b, {\bf B}) \end{align}$$ This integral can be factorized into products of three integrals over different Cartesian directions $$\begin{align} \langle G_{a} | G_{b} \rangle = S^x_{ij} S^y_{kl} S^z_{mn} \end{align}$$ Each of the integrals on the right-hand side of this equation can be obtained via the recursion relations $$\begin{align} S^{x}_{ij} &= X_{PA} S^{x}_{i-1, j} + \frac{1}{2p} \left((i-1)S^{x}_{i-2,j} + jS^{x}_{i-1,j-1}\right) \\ S^{x}_{ij} &= X_{PB} S^{x}_{i,j-1} + \frac{1}{2p} \left (iS^{x}_{i-1,j-1} + (j-1)S^{x}_{i,j-2}\right) \end{align}$$ with identical expressions for the $y$- and $z$-components of the integral. The first recursion relation can be used for any $i > 1$. The second recursion relation can be used for $i = 0$ and $j > 1$. The overlap integral corresponding to $i = j = 0$ that seeds the recursion relations is $$\begin{align} S^{x}_{00} = \sqrt{\frac{\pi}{p}}K_{ab}^{x} \end{align}$$
Now, before we try to implement these integrals, we should consider two points. First, these integrals are overlaps of primitive GTOs. As such, we will need to accumulate integrals weighted by the contraction coefficients to generate the overal overlaps between the full GTO basis fuctions. Second, we will generate the $S^{\xi}_{ij}$ workspaces for entire shell pairs, not between single basis functions, so we need to think a bit about book-keeping. As an example, for two $p$-type shells, we generate the $S^{\xi}_{ij}$ workspaces for evaluating nine different overalp integrals simultaneously. The overlap of two $p_x$-type orbitals would be $$\begin{align} \langle G_{a} | G_{b} \rangle = S^x_{11} S^y_{00} S^z_{00} \end{align}$$ whereas the overlap between $p_x$-type orbital centered at A and a $p_y$-type orbital centered at B would be $$\begin{align} \langle G_{a} | G_{b} \rangle = S^x_{10} S^y_{01} S^z_{00} \end{align}$$ We need a way to map each of these instances to a basis function with a global index. The following function will return a list of labels for a given total angular momentum $l$, ordered in a way that is consistent with the Psi4 ordering.
def get_cartesian_labels(l):
"""
generate a list of Cartesian labels corresponding to
total angular momentum, l
:param l: the total angular momentum
:return labels: the list if labels
"""
labels = []
# priority: x, then y, then z
for x in range(l, -1, -1):
for y in range(l - x, -1, -1):
z = l - y - x
labels.append((x, y, z))
return labels
Using this function, we can see that the six possible combinations of angular momentum values for $d$-type functions would be
labels = get_cartesian_labels(2)
print(labels)
[(2, 0, 0), (1, 1, 0), (1, 0, 1), (0, 2, 0), (0, 1, 1), (0, 0, 2)]
which correspond to $d_{xx}$, $d_{xy}$, $d_{xz}$, etc.
The following code loops over shell pairs and computes the overlap between the basis functions within each shell using the Obara-Saika recursion relations.
import numpy as np
# the overlap integral matrix
S = np.zeros([nbf, nbf])
# loop over shell A
for A in range(nshell):
# grab a psi4 shell object from the basis
shell_A = basis.shell(A)
# determine which atom corresponds to center A
idx_A = shell_A.ncenter
# determine the cartesisan coordinates for center A
r_A = np.array([mol.xyz(idx_A)[i] for i in range(3)])
# what is the total angular momentum for shell A?
i_A = shell_A.am
# determine the starting index for basis functions in this shell
start_idx_A = basis.shell_to_basis_function(A)
# loop over shell B
for B in range (nshell):
# grab a psi4 shell object from the basis
shell_B = basis.shell(B)
# determine which atom corresponds to center B
idx_B = shell_B.ncenter
# determine the cartesisan coordinates for center B
r_B = np.array([mol.xyz(idx_B)[j] for j in range(3)])
# what is the total angular momentum for shell B?
i_B = shell_B.am
# determine the starting index for basis functions in this shell
start_idx_B = basis.shell_to_basis_function(B)
# evaluate the distance between centers
r_AB = r_A - r_B # 9.2.14
# loop over primitives in shell A
for primitive_A in range (shell_A.nprimitive):
# exponent
exp_A = shell_A.exp(primitive_A)
# contraction coefficient
coeff_A = shell_A.coef(primitive_A)
# loop over primitives in shell B
for primitive_B in range (shell_B.nprimitive):
# exponent
exp_B = shell_B.exp(primitive_B)
# contraction coefficient
coeff_B = shell_B.coef(primitive_B)
p_AB = exp_A + exp_B # 9.2.11
mu = exp_A * exp_B / p_AB # 9.2.12
P = (exp_A * r_A + exp_B * r_B) / p_AB # 9.2.13
K_AB = np.exp(-mu * r_AB**2) # 9.2.15
# construct workspace for this shell
S_workspace = np.zeros([3, i_A+1, i_B+1]) # xyz/l/m
# seed the recursion relations with the 00 values
S_workspace[:,0,0] = (np.pi / (exp_A + exp_B))**0.5 * K_AB
# recursion ... 9.3.8, 9.3.9
r_PA = P - r_A # 9.2.14
r_PB = P - r_B # 9.2.14
for i in range (i_A+1):
for j in range (i_B+1):
if i > 0:
# 9.3.8
S_workspace[:, i, j] = r_PA * S_workspace[:, i-1, j]
if i > 1:
S_workspace[:, i, j] += (i-1) * S_workspace[:, i-2, j] / (2 * p_AB)
if j > 0:
S_workspace[:, i, j] += j * S_workspace[:, i-1, j-1] / (2 * p_AB)
elif j > 0:
# 9.3.9
S_workspace[:, i, j] = r_PB * S_workspace[:, i, j-1]
if i > 0:
S_workspace[:, i, j] += i * S_workspace[:, i-1, j-1] / (2 * p_AB)
if j > 1:
S_workspace[:, i, j] += (j-1) * S_workspace[:, i, j-2] / (2 * p_AB)
# now, we have Sx, Sy, Sz up to the desired orders
# get cartesian labels for shell A
labels_A = get_cartesian_labels(i_A)
# get cartesian labels for shell B
labels_B = get_cartesian_labels(i_B)
# loop over different angular momentum possibilities
for count_A, l_A in enumerate(labels_A):
for count_B, l_B in enumerate(labels_B):
# <a|b> = S_{ij} S_{kl} S_{mn}
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
# accumulate result with contraction coefficients
S[start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
The following Python code compares our overlap integrals to those generated by Psi4.
# get integrals from MintsHelper
mints = psi4.core.MintsHelper(wfn.basisset())
# overlap integrals
psi4_S = np.asarray(mints.ao_overlap())
print('||S - S(psi4)|| = %20.12f' % (np.linalg.norm(S - psi4_S)))
||S - S(psi4)|| = 0.000000000000
Perfect!
Multipole Moment Integrals¶
Let us now consider multipole moment integrals of the form $$\begin{align} S^{efg}_{ab} = \langle G_a | x_C^e y_C^f z_C|G_b\rangle \end{align}$$ Here, $C$ represents the origin for the multipole moments, which we will just take to be $C = [0, 0, 0]$. Like the overlap integrals, the Cartesian multipole moment integrals are separable as $$\begin{align} S^{efg}_{ab} = S_{ij}^{x,e} S_{kl}^{y,f} S_{mn}^{z,g} \end{align}$$ where $$\begin{align} S_{ij}^{x,e} = \langle G_i(x, a, A_x) | x_C^e | G_j(x, b, B_x)\rangle \end{align}$$ and we have similar expressions for the $y$- and $z$-components of the integral. Each Cartesian component of the multipole moment integrals can be obtained via the recursion relations $$\begin{align} S^{x,e}_{ij} &= X_{PA} S^{x,e}_{i-1, j} + \frac{1}{2p} \left((i-1)S^{x,e}_{i-2,j} + jS^{x,e}_{i-1,j-1} + eS^{x,e-1}_{ij} \right) \\ S^{x,e}_{ij} &= X_{PB} S^{x,e}_{i,j-1} + \frac{1}{2p} \left (iS^{x,e}_{i-1,j-1} + (j-1)S^{x,e}_{i,j-2} + eS^{x,e-1}_{ij}\right) \\ S_{ij}^{x,e+1} &= X_{PC} S_{ij}^{x,e} + \frac{1}{2p} \left ( i S_{i-1,j}^{x,e} + jS_{i,j-1}^e + eS_{ij}^{e-1}\right ) \end{align}$$ For $e=0$, we can see that the first two equations are identical to the ones we used earlier for the overlap integrals. So, for dipole integrals, we first build up the workspaces for $e=0$, followed by $e=1$. For quadrupole integrals, we then build the workspace for $e=2$.
The following code evaluates dipole and quadrupole integrals according to these recursion relations.
# dipole integrals
dipole = np.zeros([3, nbf, nbf]) # x, y, z
# quadrupole integrals
quadrupole = np.zeros([6, nbf, nbf]) # xx, xy, xz, yy, yz, zz
# origin
r_C = np.array([0, 0, 0])
# loop over shell A
for A in range(nshell):
# grab a psi4 shell object from the basis
shell_A = basis.shell(A)
# determine which atom corresponds to center A
idx_A = shell_A.ncenter
# determine the cartesisan coordinates for center A
r_A = np.array([mol.xyz(idx_A)[i] for i in range(3)])
# what is the total angular momentum for shell A?
i_A = shell_A.am
# determine the starting index for basis functions in this shell
start_idx_A = basis.shell_to_basis_function(A)
# loop over shell B
for B in range (nshell):
# grab a psi4 shell object from the basis
shell_B = basis.shell(B)
# determine which atom corresponds to center B
idx_B = shell_B.ncenter
# determine the cartesisan coordinates for center B
r_B = np.array([mol.xyz(idx_B)[j] for j in range(3)])
# what is the total angular momentum for shell A?
i_B = shell_B.am
# determine the starting index for basis functions in this shell
start_idx_B = basis.shell_to_basis_function(B)
# evaluate the distance between centers
r_AB = r_A - r_B # 9.2.14
# loop over primitives in shell A
for primitive_A in range (shell_A.nprimitive):
exp_A = shell_A.exp(primitive_A)
coeff_A = shell_A.coef(primitive_A)
# loop over primitives in shell B
for primitive_B in range (shell_B.nprimitive):
exp_B = shell_B.exp(primitive_B)
coeff_B = shell_B.coef(primitive_B)
p_AB = exp_A + exp_B # 9.2.11
mu = exp_A * exp_B / p_AB # 9.2.12
P = (exp_A * r_A + exp_B * r_B) / p_AB # 9.2.13
K_AB = np.exp(-mu * r_AB**2) # 9.2.15
# overlap workspace
S_workspace = np.zeros([3, i_A+1, i_B+1]) # xyz/l/m
S_workspace[:,0,0] = (np.pi / (exp_A + exp_B))**0.5 * K_AB
# overlap recursion ... 9.3.8, 9.3.9
r_PA = P - r_A # 9.2.14
r_PB = P - r_B # 9.2.14
r_PC = P - r_C # 9.2.14
for i in range (i_A+1):
for j in range (i_B+1):
if i > 0:
S_workspace[:, i, j] = r_PA * S_workspace[:, i-1, j]
if i > 1:
S_workspace[:, i, j] += (i-1) * S_workspace[:, i-2, j] / (2 * p_AB)
if j > 0:
S_workspace[:, i, j] += j * S_workspace[:, i-1, j-1] / (2 * p_AB)
elif j > 0:
S_workspace[:, i, j] = r_PB * S_workspace[:, i, j-1]
if i > 0:
S_workspace[:, i, j] += i * S_workspace[:, i-1, j-1] / (2 * p_AB)
if j > 1:
S_workspace[:, i, j] += (j-1) * S_workspace[:, i, j-2] / (2 * p_AB)
# S^0 has been built. build S^1 by recursion. 9.3.19
d_workspace = np.zeros([3, i_A+1, i_B+1]) # mu_xyz / xyz / l/m
for i in range (i_A+1):
for j in range (i_B+1):
# 9.3.19
d_workspace[:, i, j] = r_PC * S_workspace[:, i, j]
if i > 0:
d_workspace[:, i, j] += i * S_workspace[:, i-1, j] / (2 * p_AB)
if j > 0:
d_workspace[:, i, j] += j * S_workspace[:, i, j-1] / (2 * p_AB)
# S^1 has been built. build S^2 by recursion. 9.3.19
q_workspace = np.zeros([3, i_A+1, i_B+1]) # mu_xyz / xyz / l/m
for i in range (i_A+1):
for j in range (i_B+1):
# 9.3.19
q_workspace[:, i, j] = r_PC * d_workspace[:, i, j]
q_workspace[:, i, j] += S_workspace[:, i, j] / (2 * p_AB)
if i > 0:
q_workspace[:, i, j] += i * d_workspace[:, i-1, j] / (2 * p_AB)
if j > 0:
q_workspace[:, i, j] += j * d_workspace[:, i, j-1] / (2 * p_AB)
# we have dip_x, dip_y, dip_z up to desired orders
labels_A = get_cartesian_labels(i_A)
labels_B = get_cartesian_labels(i_B)
for count_A, l_A in enumerate(labels_A):
for count_B, l_B in enumerate(labels_B):
# x
prod = d_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
dipole[0, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# y
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= d_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
dipole[1, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# z
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= d_workspace[2, l_A[2], l_B[2]]
dipole[2, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# we also have q_xx, q_yy, q_zz up to desired orders
for count_A, l_A in enumerate(labels_A):
for count_B, l_B in enumerate(labels_B):
# xx
prod = q_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
quadrupole[0, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# xy
prod = d_workspace[0, l_A[0], l_B[0]]
prod *= d_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
quadrupole[1, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# xz
prod = d_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= d_workspace[2, l_A[2], l_B[2]]
quadrupole[2, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# yy
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= q_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
quadrupole[3, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# yz
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= d_workspace[1, l_A[1], l_B[1]]
prod *= d_workspace[2, l_A[2], l_B[2]]
quadrupole[4, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# zz
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= q_workspace[2, l_A[2], l_B[2]]
quadrupole[5, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# don't forget the factor of -1 from the electron charge
dipole *= -1
quadrupole *= -1
Now, we can compare our multipole moment integrals to those generated by Psi4.
dir = ['x', 'y', 'z']
# dipole integrals
psi4_dipole = np.asarray(mints.ao_dipole())
for i in range (3):
print('||mu(%s) - mu(%s,psi4)|| = %20.12f' % (dir[i], dir[i], np.linalg.norm(dipole - psi4_dipole)))
print('')
# quadrupole integrals
psi4_quadrupole = np.asarray(mints.ao_quadrupole())
count = 0
for i in range (3):
for j in range (i, 3):
print('||q(%s%s) - q(%s%s,psi4)|| = %20.12f' % (dir[i], dir[j], dir[i], dir[j], np.linalg.norm(quadrupole[count] - psi4_quadrupole[count])))
count += 1
||mu(x) - mu(x,psi4)|| = 0.000000000000 ||mu(y) - mu(y,psi4)|| = 0.000000000000 ||mu(z) - mu(z,psi4)|| = 0.000000000000 ||q(xx) - q(xx,psi4)|| = 0.000000000000 ||q(xy) - q(xy,psi4)|| = 0.000000000000 ||q(xz) - q(xz,psi4)|| = 0.000000000000 ||q(yy) - q(yy,psi4)|| = 0.000000000000 ||q(yz) - q(yz,psi4)|| = 0.000000000000 ||q(zz) - q(zz,psi4)|| = 0.000000000000
Linear Momentum and Kinetic Energy Integrals¶
Integrals involving the linear momentum operator $$\begin{align} \hat{p}_x = -i\frac{\partial}{\partial x} \end{align}$$ etc. and the kinetic energy operator $$\begin{align} \hat{T} = -\frac{1}{2}\nabla^2 \end{align}$$ can be evaluated using Obara-Saika recursion relations for differential operators. As above, integrals over differential operators of the form $$\begin{align} D^{efg}_{ab} = \langle G_a | \frac{\partial^e}{\partial x^e}\frac{\partial^f}{\partial y^f}\frac{\partial^g}{\partial z^g} | G_b \rangle \end{align}$$ can be factorized along each Cartesian direction as $$\begin{align} D^{efg}_{ab} = D^{x,e}_{ij} D^{y,f}_{kl} D^{z,g}_{mn} \end{align}$$ Each component can be obtained from the recursion relations $$\begin{align} D^{x,e}_{ij} &= X_{PA} D^{x,e}_{i-1, j} + \frac{1}{2p} \left((i-1)D^{x,e}_{i-2,j} + jD^{x,e}_{i-1,j-1} - 2beD^{x,e-1}_{i-1,j} \right) \\ D^{x,e}_{ij} &= X_{PB} D^{x,e}_{i,j-1} + \frac{1}{2p} \left (iD^{x,e}_{i-1,j-1} + (j-1)D^{x,e}_{i,j-2} + 2aeD^{x,e-1}_{i,j-1}\right) \\ D_{ij}^{x,e+1} &= 2a D_{i+1,j}^{x,e} - i D_{i-1, j}^{x,e} \end{align}$$ As we saw with the multipole moment integrals, the case of $e=0$ is special; the recursion relations reduce to those for the overlap integrals. So, we can start by constructing the workspace for $e=0$, and then we can build the integrals for $e=1$ using the third recursion relation. Note that each application of a differential operator ($e+1$) will require knowledge of one higher unit of angular momentum at the lower level ($e$). For example, if second derivative integrals are desired, the overlap integral workspace should be constructed to two units of angular momentum higher than the angular momentum of the primitive GTOs in that shell.
The following code evaluates momentum and kinetic energy integrals according to these recursion relations.
# kinetic energy
T = np.zeros([nbf, nbf])
# nabla (the derivative part of the momentum operator)
nabla = np.zeros([3, nbf, nbf])
# loop over shell A
for A in range(nshell):
# grab a psi4 shell object from the basis
shell_A = basis.shell(A)
# determine which atom corresponds to center A
idx_A = shell_A.ncenter
# determine the cartesisan coordinates for center A
r_A = np.array([mol.xyz(idx_A)[i] for i in range(3)])
# what is the total angular momentum for shell A?
i_A = shell_A.am
# determine the starting index for basis functions in this shell
start_idx_A = basis.shell_to_basis_function(A)
# loop over shell B
for B in range (nshell):
# grab a psi4 shell object from the basis
shell_B = basis.shell(B)
# determine which atom corresponds to center B
idx_B = shell_B.ncenter
# determine the cartesisan coordinates for center B
r_B = np.array([mol.xyz(idx_B)[j] for j in range(3)])
# what is the total angular momentum for shell A?
i_B = shell_B.am
# determine the starting index for basis functions in this shell
start_idx_B = basis.shell_to_basis_function(B)
# evaluate the distance between centers
r_AB = r_A - r_B # 9.2.14
# loop over primitives for shell A
for primitive_A in range (shell_A.nprimitive):
exp_A = shell_A.exp(primitive_A)
coeff_A = shell_A.coef(primitive_A)
# loop over primitives for shell B
for primitive_B in range (shell_B.nprimitive):
exp_B = shell_B.exp(primitive_B)
coeff_B = shell_B.coef(primitive_B)
p_AB = exp_A + exp_B # 9.2.11
mu = exp_A * exp_B / p_AB # 9.2.12
P = (exp_A * r_A + exp_B * r_B) / p_AB # 9.2.13
K_AB = np.exp(-mu * r_AB**2) # 9.2.15
S_workspace = np.zeros([3, i_A+3, i_B+3]) # xyz/l/m
S_workspace[:,0,0] = (np.pi / (exp_A + exp_B))**0.5 * K_AB
# recursion ... 9.3.8, 9.3.9
r_PA = P - r_A # 9.2.14
r_PB = P - r_B # 9.2.14
for i in range (i_A+3):
for j in range (i_B+3):
if i > 0:
# 9.3.8 ... with i+1 -> i, i -> i-1
S_workspace[:, i, j] = r_PA * S_workspace[:, i-1, j]
if i > 1:
S_workspace[:, i, j] += (i-1) * S_workspace[:, i-2, j] / (2 * p_AB)
if j > 0:
S_workspace[:, i, j] += j * S_workspace[:, i-1, j-1] / (2 * p_AB)
if j > 0:
# 9.3.9 ... with j+1 -> j, j -> j-1
S_workspace[:, i, j] = r_PB * S_workspace[:, i, j-1]
if i > 0:
S_workspace[:, i, j] += i * S_workspace[:, i-1, j-1] / (2 * p_AB)
if j > 1:
S_workspace[:, i, j] += (j-1) * S_workspace[:, i, j-2] / (2 * p_AB)
# now, S is built. we can build nabla by recursion ... 9.3.30
nabla_workspace = np.zeros([3, i_A+2, i_B+2]) # p_xyz / l/m
for i in range (i_A+2):
for j in range (i_B+2):
# 9.3.30
nabla_workspace[:, i, j] = 2 * exp_A * S_workspace[:, i+1, j]
if i > 0:
nabla_workspace[:, i, j] -= i * S_workspace[:, i-1, j]
# we can build T by recursion ... 9.3.31
T_workspace = np.zeros([3, i_A+1, i_B+1]) # T_xyz / xyz / l/m
for i in range (i_A+1):
for j in range (i_B+1):
# 9.3.31
T_workspace[:, i, j] = 4 * exp_A**2 * S_workspace[:, i+2, j]
T_workspace[:, i, j] -= 2 * exp_A * (2 * i + 1) * S_workspace[:, i, j]
if i > 1:
T_workspace[:, i, j] += i * (i - 1) * S_workspace[:, i-2, j]
# now ... we have nabla_x, nabla_y, nabla_z up to desired orders
labels_A = get_cartesian_labels(i_A)
labels_B = get_cartesian_labels(i_B)
for count_A, l_A in enumerate(labels_A):
for count_B, l_B in enumerate(labels_B):
prod = nabla_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
nabla[0, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= nabla_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
nabla[1, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= nabla_workspace[2, l_A[2], l_B[2]]
nabla[2, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# now ... we have Tx, Ty, Tz up to desired orders
for count_A, l_A in enumerate(labels_A):
for count_B, l_B in enumerate(labels_B):
prod = T_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
T[start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= T_workspace[1, l_A[1], l_B[1]]
prod *= S_workspace[2, l_A[2], l_B[2]]
T[start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
prod = S_workspace[0, l_A[0], l_B[0]]
prod *= S_workspace[1, l_A[1], l_B[1]]
prod *= T_workspace[2, l_A[2], l_B[2]]
T[start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * prod
# don't forget the factor of -1/2 on the kinetic energy integrals
T *= -0.5
Now, we can compare our derivative integrals to those generated by Psi4.
# nabla integrals
psi4_nabla = np.asarray(mints.ao_nabla())
for i in range (3):
print('||nabla(%s) - nabla(%s,psi4)|| = %20.12f' % (dir[i], dir[i], np.linalg.norm(nabla - psi4_nabla)))
print('')
# kinetic energy integrals
psi4_T = np.asarray(mints.ao_kinetic())
print('||T - T(psi4)|| = %20.12f' % (np.linalg.norm(T - psi4_T)))
||nabla(x) - nabla(x,psi4)|| = 0.000000000000 ||nabla(y) - nabla(y,psi4)|| = 0.000000000000 ||nabla(z) - nabla(z,psi4)|| = 0.000000000000 ||T - T(psi4)|| = 0.000000000000
Angular Momentum Integrals¶
Lastly, we can use the Obara-Saika schemes for constructing dipole integrals and integrals over differential operators to generate integrals over the angular momentum operator $$\begin{align} \bf{L} = {\bf r} \times {\bf p} \end{align}$$ which has Cartesian components $$\begin{align} L_x &= -i \left ( y \frac{\partial}{\partial z} - z \frac{\partial}{\partial y} \right ) \\ L_y &= -i \left ( z \frac{\partial}{\partial x} - x \frac{\partial}{\partial z} \right ) \\ L_z &= -i \left ( x \frac{\partial}{\partial y} - y \frac{\partial}{\partial x} \right ) \\ \end{align}$$ The integral $$\begin{align} i\langle G_a | L_z | G_b\rangle = \langle G_a | x \frac{\partial}{\partial y} - y \frac{\partial}{\partial x} | G_b \rangle \end{align}$$ factorizes into a product of an overlap integral in the $z$-direction and multipole moment and kinetic energy integrals in the $x$- and $y$-directions: $$\begin{align} \langle G_a | x \frac{\partial}{\partial y} - y \frac{\partial}{\partial x} | G_b \rangle = (S^{x,1}_{ij} D^{y,1}_{kl} - D^{x,1}_{ij} S^{y,1}_{kl}) S^z_{mn} \end{align}$$ The following code implements the angular momentum integrals in this way.
# angular momentum integrals
am = np.zeros([3, nbf, nbf])
# loop over shell A
for A in range(nshell):
# grab a psi4 shell object from the basis
shell_A = basis.shell(A)
# determine which atom corresponds to center A
idx_A = shell_A.ncenter
# determine the cartesisan coordinates for center A
r_A = np.array([mol.xyz(idx_A)[i] for i in range(3)])
# what is the total angular momentum for shell A?
i_A = shell_A.am
# determine the starting index for basis functions in this shell
start_idx_A = basis.shell_to_basis_function(A)
# loop over shell B
for B in range (nshell):
# grab a psi4 shell object from the basis
shell_B = basis.shell(B)
# determine which atom corresponds to center B
idx_B = shell_B.ncenter
# determine the cartesisan coordinates for center B
r_B = np.array([mol.xyz(idx_B)[j] for j in range(3)])
# what is the total angular momentum for shell A?
i_B = shell_B.am
# determine the starting index for basis functions in this shell
start_idx_B = basis.shell_to_basis_function(B)
# evaluate the distance between centers
r_AB = r_A - r_B # 9.2.14
# loop over primitives for shell A
for primitive_A in range (shell_A.nprimitive):
exp_A = shell_A.exp(primitive_A)
coeff_A = shell_A.coef(primitive_A)
# loop over primitives for shell B
for primitive_B in range (shell_B.nprimitive):
exp_B = shell_B.exp(primitive_B)
coeff_B = shell_B.coef(primitive_B)
p_AB = exp_A + exp_B # 9.2.11
mu = exp_A * exp_B / p_AB # 9.2.12
P = (exp_A * r_A + exp_B * r_B) / p_AB # 9.2.13
K_AB = np.exp(-mu * r_AB**2) # 9.2.15
S_workspace = np.zeros([3, i_A+3, i_B+3]) # xyz/l/m
S_workspace[:,0,0] = (np.pi / (exp_A + exp_B))**0.5 * K_AB
# recursion ... 9.3.8, 9.3.9
r_PA = P - r_A # 9.2.14
r_PB = P - r_B # 9.2.14
r_PC = P - r_C # 9.2.14
for i in range (i_A+2):
for j in range (i_B+2):
if i > 0:
S_workspace[:, i, j] = r_PA * S_workspace[:, i-1, j]
if i > 1:
S_workspace[:, i, j] += (i-1) * S_workspace[:, i-2, j] / (2 * p_AB)
if j > 0:
S_workspace[:, i, j] += j * S_workspace[:, i-1, j-1] / (2 * p_AB)
if j > 0:
S_workspace[:, i, j] = r_PB * S_workspace[:, i, j-1]
if i > 0:
S_workspace[:, i, j] += i * S_workspace[:, i-1, j-1] / (2 * p_AB)
if j > 1:
S_workspace[:, i, j] += (j-1) * S_workspace[:, i, j-2] / (2 * p_AB)
# S has been built. we can build nabla by recursion ... 9.3.30
nabla_workspace = np.zeros([3, i_A+1, i_B+1]) # p_xyz / l/m
for i in range (i_A+1):
for j in range (i_B+1):
# 9.3.30
nabla_workspace[:, i, j] = 2 * exp_A * S_workspace[:, i+1, j]
if i > 0:
nabla_workspace[:, i, j] -= i * S_workspace[:, i-1, j]
# S^0 has been built. now build S^1 by recursion ... 9.3.19
d_workspace = np.zeros([3, i_A+1, i_B+1]) # mu_xyz / xyz / l/m
for i in range (i_A+1):
for j in range (i_B+1):
# 9.3.19
d_workspace[:, i, j] = r_PC * S_workspace[:, i, j]
if i > 0:
d_workspace[:, i, j] += i * S_workspace[:, i-1, j] / (2 * p_AB)
if j > 0:
d_workspace[:, i, j] += j * S_workspace[:, i, j-1] / (2 * p_AB)
# we have nabla_x, nabla_y, nabla_z, mu_x, etc. up to desired orders
labels_A = get_cartesian_labels(i_A)
labels_B = get_cartesian_labels(i_B)
for count_A, l_A in enumerate(labels_A):
for count_B, l_B in enumerate(labels_B):
# y d/dz - z d/dy
prod1 = S_workspace[0, l_A[0], l_B[0]]
prod1 *= d_workspace[1, l_A[1], l_B[1]]
prod1 *= nabla_workspace[2, l_A[2], l_B[2]]
prod2 = S_workspace[0, l_A[0], l_B[0]]
prod2 *= nabla_workspace[1, l_A[1], l_B[1]]
prod2 *= d_workspace[2, l_A[2], l_B[2]]
am[0, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * (prod1 - prod2)
# z d/dx - x d/dz
prod1 = nabla_workspace[0, l_A[0], l_B[0]]
prod1 *= S_workspace[1, l_A[1], l_B[1]]
prod1 *= d_workspace[2, l_A[2], l_B[2]]
prod2 = d_workspace[0, l_A[0], l_B[0]]
prod2 *= S_workspace[1, l_A[1], l_B[1]]
prod2 *= nabla_workspace[2, l_A[2], l_B[2]]
am[1, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * (prod1 - prod2)
# x d/dy - y d/dx
prod1 = d_workspace[0, l_A[0], l_B[0]]
prod1 *= nabla_workspace[1, l_A[1], l_B[1]]
prod1 *= S_workspace[2, l_A[2], l_B[2]]
prod2 = nabla_workspace[0, l_A[0], l_B[0]]
prod2 *= d_workspace[1, l_A[1], l_B[1]]
prod2 *= S_workspace[2, l_A[2], l_B[2]]
am[2, start_idx_A + count_A, start_idx_B + count_B] += coeff_A * coeff_B * (prod1 - prod2)
# don't forget the factor of -1
am *= -1
Now, we can compare our angular momentum integrals to those generated by Psi4.
# check angular momentum integrals
psi4_am = np.asarray(mints.ao_angular_momentum())
for i in range (3):
print('||am(%s) - am(%s,psi4)|| = %20.12f' % (dir[i], dir[i], np.linalg.norm(am - psi4_am)))
print('')
||am(x) - am(x,psi4)|| = 0.000000000000 ||am(y) - am(y,psi4)|| = 0.000000000000 ||am(z) - am(z,psi4)|| = 0.000000000000