4.1. Hartree-Fock Procedure for Approximate Quantum Chemistry#
import psi4
import numpy as np
from scipy import linalg as splinalg
import sys
sys.path.append("..")
from helpers import *
We learned that we can only solve the Schroedinger equation exactly for one-electron systems.
To solve many-electron systems we need to use Hartree-Fock molecular orbital theory. Just as a test case, let’s use Psi4 to compute the Hartree-Fock wavefunction and energy for the Hydrogen atom:
In this lab activity, you will build and diagonalize the Fock matrix to determine the MO coefficients and energies for a molecule. We will be using the functions of the Psi4 quantum chemistry software package to compute the integrals we need. The following notebook will lead you through setting up your molecule, establishing the basis set, and forming and diagonalizing the Fock matrix. Be sure to run each cell as your proceed through the notebook.
4.1.1. The Hartree-Fock procedure#
The Schrödinger equation has the structure of an eigenvalue equation
In Hartree-Fock theory, this is reexpresed in terms of the Fock matrix, \(F\), a matrix of wavefunction amplitudes for each MO, \(C\), and the overlap matrix, \(S\),
The Fock matrix for a closed-shell system is
where \(H\) is the one electron “core” Hamiltonian, \(J\) is the Coulomb integral matrix, and \(K\) is the exchange integral matrix. The definitions of \(J\) and \(K\) depend on the coefficients, \(C\). We see here the central premise of SCF: To get the Fock matrix, we need the coefficient matrix, but to compute the coefficient matrix we need the Fock matrix. When \(S\) is not equal to the identity matrix (i.e. the basis is not orthonormal), then this is a pseudo-eigenvalue problem and is even harder to solve. This is our task.
We will
introduce the overlap matrix
learn how to build an orthogonalization matrix
learn how to calculate the density
learn how to calculate the Coulomb and Exchange integral matrices
learn how to diagonalize the Fock matrix
build an iterative procedure to converge HF energy
It can be helpful to look back at exercise1
4.1.2. Orthogonalizing the AO basis set#
Use the two orthonormal basis vector that we defined in Exercise 1
# define a first basis vector and a second, orthogonal vector (you can use your solutions from Ex 1)
phi1 = np.array([]) #insert here
phi2 = np.array([]) #insert here
4.1.3. IV. The overlap integrals#
For a set of basis functions, \(\phi_i(\tau)\), where \(\tau\) is a shorthand for all the coordinates of all the particles, we can calculate the overlap integrals between the basis functions in the following way
Exercise 1
Define \(S_{ij}\) using Dirac notation.
Exercise 2
For an orthonormal basis, what does the overlap integral array, S, look like?
Calculate the terms S_ij using the basis vectors phi1 and phi2.
# calculate the overlap (inner product) of the vectors
S_11 = None #<your formulas here>
S_12 = None
S_21 = None
S_22 = None
print('The ij elements of S:')
print(S_11,S_12)
print(S_21,S_22)
4.1.4. Constructing the overlap matrix#
These overlap integrals, \(S_{ij}\), can be interpreted as the elements on the \(i\)-th row and \(j\)-th column of a matrix, \(S\). Let’s propose a matrix, \(S\), made of the overlap integrals \(S_{ij}\). We can build \(S\) systematically in the following way. First, make a matrix, \(B\), composed of our basis vectors as columns,
We will use the symbol \(\dagger\) to indicate the complex conjugate of the transpose of a matrix. So
# construct the overlap matrix from matrix of basis vectors
vector_length = phi1.size #length of the vector space
phi1_column = phi1.reshape(vector_length,1) #this makes phi a column vector
phi2_column = phi2.reshape(vector_length,1)
# put together (concatenate) the vectors into the matrix B
B = np.concatenate((phi1_column,phi2_column),axis=1)
print(F'The matrix B:\n{B}')
B_dagger = B.conj().T
print(F'The matrix B^\\dagger:\n{B_dagger}')
Now, multiplying the rows of \(B^\dagger\) by the columns of \(B\) (normal matrix multiplication) produces a matrix of the overlap integrals in the correct locations. we defined to be the matrix \(S\).
Exercise 3
Use B and B_dagger and the matrix rules above to calculate the matrix S.
# calculate S from matrix of basis vectors
S = None #<your formula here>
print(F'The matrix of eigenvectors in columns B =\n {B} \n\nand S = B^† B =\n{S}\n')
4.1.5. Einstein implicit summation notation#
Matrix multiplication is defined through
Note that there is a repeated index, \(i\), in the summation. In implicit summation notation, Einstein notation, we do not write the \(\sum\) and treat the summation as understood.
Using implicit summation for the case at hand, \(B^\dagger B\) gives
where \(B^*\) is the complex conjugate of \(B\) (no transpose). Note that the two sets of indices, \((i,p)\) and \((i,q)\), in the input matrices become one set, \((p,q)\), in the product.
Note
There is a convenient function in numpy called einsum(), which is one of the crown jewels of the numpy library. In short, einsum lets you perform various combinations of multiplying, summing, and transposing matrices very efficiently. This is a good tutorial about einsum if you want to learn more. In this exercise, it suffices to look carefully at the indices of the matrices and use the Einstein Summation convention.)
So in order to calculate the sum
use
S=numpy.einsum('ip,iq->pq',B.conj(),B)
We will use einsum() in several places.
Let’s try this with a simple basis set of two (perhaps) orthonormal vectors.
Exercise 4
Describe how the notation of the np.einsum command correlates to the implicit summation formula written above.
Exercise 5
Use the function np.einsum() to calculate the matrix S, and confirm that your answer is the same as above.
# calculate S from Einstein sum
S = None #<your formula here>
print(F'S from Einstein notation:\n{S}')
Exercise 6
Propose a different orthonormal basis, modify phi1 and phi2, and verify that S still has the same form. There are infinitely many choices. It isn’t complex… or is it?!
4.1.6. Gaussian atomic orbital basis set#
H\(_2\)O is a small but interesting molecule to use in our exploration.
Exercise 7
How many electrons are there in total in H\(_2\)O? How many occupied molecular orbitals would you expect?
Before we can begin to implement the HF procedure, we need to specify the molecule and basis set that we will be using. We will also set the memory usage for our calcluation and the output file name.
# ==> Set Basic Psi4 Options <==
# Memory specification
psi4.set_memory('500 MB')
numpy_memory = 2 # No NumPy array can exceed 2 MB in size
# set output file
psi4.core.set_output_file('output.dat', False)
# specify the basis
# basis = 'cc-pvdz'
basis = 'sto-3g'
# Set computation options
psi4.set_options({'basis': basis,
'scf_type': 'pk',
'e_convergence': 1e-8})
# ==> Define Molecule <==
# Define our model of water
mol = psi4.geometry("""
symmetry c1
O
H 1 1.1
H 1 1.1 2 104
""")
# compute energy
SCF_E_psi = psi4.energy('scf')
psi4.core.clean()
print(f"The Hartree-Fock ground state energy of the water is: {SCF_E_psi} Eh")
Memory set to 476.837 MiB by Python driver.
The Hartree-Fock ground state energy of the water is: -74.94207989868094 Eh
Next, we need to build the wavefunction from the basis functions. We store the wavefunction in a variable called wfn. We use the function nalpha() provided by the wavefunction object we created above, wfn, to determine the number of orbitals with spin alpha, which will be doubly occupied orbitals for close shelled systems. We save this answer as a variable called ndocc (number of doubly occupied orbitals).
Execute the code below and confirm that the number of doubly occupied orbitals matches your expectation for the molecule you chose.
# ==> Compute static 1e- and 2e- quantities with Psi4 <==
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('basis'))
# number of spin alpha orbitals (doubly occupied for closed-shell systems)
ndocc = wfn.nalpha()
nbf = wfn.basisset().nbf()
print(F'Number of occupied orbitals: {ndocc}')
print(F'Number of basis functions: {nbf}')
Next we will examine the atomic orbital basis set. To do this, we have to set up a data structure, called a class, to calculate the molecular integrals. (Psi4 will do the nasty calculus for us.) We will call this data structure mints (Molecular INTegralS). We use the function ao_overlap to calculate the overlap integrals between all the AO basis functions. We cast the result to a numpy array called S.
# Construct a molecular integrals object
mints = psi4.core.MintsHelper(wfn.basisset())
# Overlap matrix as a psi4 Matrix object
S_matrix = mints.ao_overlap()
# Overlap matrix converted into an ndarray
S = np.asarray(S_matrix)
print(F'Shape of S is {S.shape}')
Exercise 8
Explain the shape (number of rows and columns) of S in terms of the AO basis set we chose.
Examine the contents of S.
print(S) #the full matrix may be somewhat hard to read based on the basis set
# Look at the first few elements
peak(S)
Exercise 9
Based on your observations of S in the AO basis, answer the following questions
What do the diagonal elements of
Sindicate?What do the off-diagonal elements of
Sindicate?Does the Gaussian atomic orbital basis set form an orthonormal basis?
We can perform this test programmatically as well, with a few python tricks. Construct an array of the same size as the overlap array (S) that has 1’s along the diagonal and 0’s everywhere else. Then compare that array to the S array to doubly verify whether or not you think the AO basis is orthonormal.
# use some functions from helpers library, you can test your overlap arrays here
isBasisOrthonormal(S)
Exercise 10
Does the result of your extra evaluation agree with what you determined previously?
4.1.7. An orthogonalization matrix#
Recall that if we had used the hydrogen atom wavefunctions as our basis set, the AO wavefunctions would all be orthonormal. Since we used a basis set of Gaussian wavefuctions, this may not be the case. We will now introduce some tools to fix it!
Since our AO basis set was not orthonormal, we seek to construct an orthogonalization matrix, \(A\), such that \({\bf A}^{\dagger}{\bf SA} = {\bf 1}\).
Motivation: If \({\bf A}\) and \({\bf S}\) were real numbers \(a\) and \(s\) (not matrices), this would be simple to solve. First, \(a^\dagger=a\) because a real number is the same as its hermitian transpose. By simple algebra we can solve for a,
In linear algebra (with matrices instead of numbers) this is more complicated, but numpy and the mints class can take care of the details for us! Leaving out the details, we will calculate
Exercise 11
Use the function np.linalg.inv() to calculate the inverse of S, and the function splinalg.sqrtm() to take its (matrix) square root. Execute the code below and examine the matrix A.
# ==> Construct AO orthogonalization matrix A <==
# inverse of S using np.linalg.inv
# matrix square root of the inverse of S using splinalg.sqrtm
A = None
peak(A)
Exercise 12
What do you observe about the elements of A? Is the matrix real or complex? Is the matrix symmetric or not?
Our basis set \(B\) is not orthonormal, so we want to take linear combinations of its columns to make a new basis set, \(B'\), that is orthonormal. We define a new matrix, \(A\), in terms of that transformation,
The new overlap matrix will be
The matrix \(A\) makes the proper linear combination of the columns of \(B\) and \(A^\dagger\) makes the linear combinations of the rows of \(B^\dagger\). This is a very common structure of transformation matrices. Because \(S\) is real and symmetric, \(A\) is also real and symmetric, so \(A^\dagger=A\). The transformation becomes simply
Exercise 13
Use the orthogonalization matrix A to transform the overlap matrix, S. Check the transformed overlap matrix, S_p, to make sure it represents an orthonormal basis.
# Transform S with A and assign the result to S_p
S_p = None #<your equation here>
isBasisOrthonormal(S_p)
Exercise 14
The product A S A does not take the complex conjugate transpose of A. What conditions (properties of A) make that ok?
4.1.8. The Fock Matrix Transformed to an Orthonormal Basis#
We can now return to Eq. (4.1) using our orthogonalization matrix \(A\). A common linear algebra trick is to “insert one.” In this case, the matrix \(A\) times its inverse is, by definition the identity matrix, \(AA^{-1}={\bf1}\). We can put that factor of one anywhere in an equation that is useful to us, and, then, typically we regroup terms in a way we want. In this case, we insert one between \(FC\) and \(SC\).
Multiplying on the left by \(A\) then gives
We can recognize the transformation \(S'=ASA\) on the right hand side and similarly define \(F'=AFA\) on the left hand side. Lastly, we define a transformed coefficient matrix, \(C'=A^{-1}C\). Our transformed Fock equation reads
The last line follows because \(S'={\bf1}\) in our new basis. We now have an eigenvalue problem that we can solve by matrix diagonalization.
In the expression
the matrix \(A^{-1}\) transforms the coefficients, \(C\), into the orthogonalized basis set. We will also need a way to transform those coefficients, \(C'\), back to the original AO basis.
Exercise 15
Based on the definition of \(C'\), propose a definition of \(C\) in terms of \(A\) and \(C'\). Justify your equation.
4.1.9. Initial guess for the Fock Matrix is the one electron Hamiltonian#
To get the Fock matrix, we need the coefficient matrix, but to compute the coefficient matrix we need the Fock matrix. So we start with a guess for the Fock matrix, which is the core Hamiltonian matrix.
# Build core Hamiltonian
T = np.asarray(mints.ao_kinetic())
V = np.asarray(mints.ao_potential())
H = T + V
Exercise 16
In the cell below, use the core Hamiltonian matrix as your initial guess for the Fock matrix. Transform it with the same A matrix you used above. To calculate the eigenvalues, vals, and eigenvectors, vecs, of matrix M using vals, vecs = np.linalg.eigh(M).
# Guess for the Fock matrix
F = None # Replace
# Transformed Fock matrix
F_p = None # Replace
# Diagonalize F_p for eigenvalues & eigenvectors with NumPy
e, C_p = None # Replace
Exercise 17
Display, i.e., print, the coefficent matrix and confirm it the correct size
print(C_p)
Now that we have the coefficents in the transformed basis, we need to go back and get the coefficients in the original AO basis.
Exercise 18
Use A and the formula you proposed previously to transform the coefficient matrix back to the AO basis. Confirm that the resulting matrix appears reasonable, i.e., similar size and magnitude
# Transform the coefficient matrix back into AO basis
C = A.dot(C_p)
print(C)
4.1.10. The Density Matrix#
Recall, the Fock matrix is
where \(H\) is the one electron “core” Hamiltonian, \(J\) is the Coulomb integral matrix, and \(K\) is the exchange integral matrix. The HF energy can be expressed in explicit terms of one and two electron integrals
Expanding the orbitals in terms of basis functions, we find
First, look at the coefficients in (4.5). They come in two pairs of complex-conjugates, \(c^*_{pi}c_{qi}\) and \(c^*_{ri}c_{si}\). The diagonal terms, when \(p=q\) for example, are the probability of some basis function \(p\) contributing to the MO \(i\). We will sum each term over the occupied orbitals, \(i\), to form the “density matrix”
We are going to construct the density matrix from the occupied orbitals. To get a matrix of just the occupied orbitals, use the coefficient matrix in the original AO basis, and take a slice to include all rows and just the columns that represent the occupied orbitals.
# Grab occupied orbitals (recall: ndocc is the number of doubly occupied orbitals we found earlier)
C_occ = C[:, :ndocc]
print(F'The shape of C_occ is {C_occ.shape}')
Exercise 19
Build the density matrix, D, from the occupied orbitals, C_occ, using the function np.einsum(). Hint Look at (4.6)
# Build density matrix from occupied orbitals
D = np.array([]) # replace this command with np.einsum() as described above
print(F'The shape of D is {D.shape}')
4.1.11. Coulomb and Exchange Integrals and the SCF Energy#
The integral on the right of (4.7) is super important. It has four indicies, \(p,q,r,s\), so formally it is a tensor. It accounts for the repulsion between pairs of electrons, so it is called the electron repulsion integral tensor, \(I\),
First, we can build the electron-repulsion integral (ERI) tensor, which calculates the electron repulsion between the atomic orbital wavefunctions, and the core Hamiltonian. Mints does all the work for us!
# Build electron repulsion integral (ERI) Tensor
I = np.asarray(mints.ao_eri())
(4.5) can be expressed in terms of \(I\) and \(D\) as
The term \(\sum_{rs}D_{rs}I_{pqrs}\) is the effective repulsion felt by one electron due to the other electrons in the system. This term is the Coulomb integral matrix
Exercise 20
Define J in terms of the density matrix, D, and the electron repulsion integral tensor, I, using np.einsum().
#Define J
J = None # use np.einsum here
Similarly,
corrects the repulsion due to electrons “avoiding each other” due to their Fermionic (antisymmetric w.r.t. exchange) character. This term is the exchange integral matrix
Exercise 21
Define K in terms of the density matrix, D, and the electron repulsion integral tensor, I, using einsum().
#Define K
K = None # use np.einsum here
Exercise 22
Define F in terms of H, J, and K. (Recall The Hartree-Fock procedure)
#Define F as a function of H, J, and K
F = None # insert code here
A more convenient form of the SCF energy is in terms of a sum over the AO basis functions
Exercise 23
Calculate the SCF energy based on H, F, and D using np.einsum().
E_nuc = mol.nuclear_repulsion_energy()
SCF_E = None # insert your code here, don't forget to use E_nuc
print(F'Energy = {SCF_E:.8f}')
Exercise 24
Based on the result of the calculation in basisset, is this a reasonable answer?
Exercise 25
Describe a procedure (i.e. identify the steps) that will update coefficients and compute a new density matrix based on the updated values of the Fock matrix.
Exercise 26
Using the procedure proposed above, calculate the updated coefficients
# insert remaining steps here, you should end with a new C
# transform F to the orthonormal basis
F_p =
# find e and C_p
e, C_p = n
# find C
C =
# Select the occupied MOs
C_occ = C[:, :ndocc]
# Update the density matrix
D = np.einsum('pi,qi->pq', C_occ, C_occ)
You have just completed one cycle of the SCF calculation!
Now we will use the density matrix to build the Fock matrix. The code block below sets up a skeleton of the Hartree-Fock procedure. The basic steps are:
Calculate the Fock Matrix based on the density matrix previously defined from a one electron hamiltonian
Calculate the energy from the Fock matrix.
Check and see if the energy has converged by comparing the current energy to the previous energy and seeing if it is within the convergence threshold.
If the energy has not converged, transform the Fock matrix, and diagonalize the transformed Fock matrix to get the energy and MO coefficients. Then transform back to the original AO basis, pull the occupied orbitals, and reconstruct the density matrix.
# ==> Nuclear Repulsion Energy <==
E_nuc = mol.nuclear_repulsion_energy()
# ==> SCF Iterations <==
# Pre-iteration energy declarations
SCF_E = 0.0
E_old = 0.0
# ==> Set default program options <==
# We continue recalculating the energy until it converges to the level we specify.
# The varible `E_conv` is where we set this level of convergence. We also set a
# maximum number of iterations so that if our calculation does not converge, it
# eventually stops and lets us know that it did not converge.
# Maximum SCF iterations
MAXITER = 40
# Energy convergence criterion
E_conv = 1.0e-6
print('==> Starting SCF Iterations <==\n')
# Begin Iterations
for scf_iter in range(1, MAXITER + 1):
# Build Fock matrix (meaning define J, K with density matrix and then define F with J and k)
# <your code here>
# Compute SCF energy (don't forget E_nuc!)
SCF_E = #<your formula here>
print(F'SCF Iteration {scf_iter}: Energy = {SCF_E:.8f} dE = {SCF_E - E_old:.8f}')
# Check to see if the energy is converged. If it is break out of the loop.
# If it is not, set the current energy E_old
if (abs(SCF_E - E_old) < E_conv):
break
E_old = SCF_E
# Compute new coefficient & density matrices (Exercise 26)
# <your code here>
# MAXITER exceeded?
if (scf_iter == MAXITER):
psi4.core.clean()
raise Exception("Maximum number of SCF iterations exceeded.")
# Post iterations
print('\nSCF converged.')
print(F'Final RHF Energy: {SCF_E:.6f} [Eh]')
Compare your results to Psi4 by computing the energy using psi4.energy() in the cell below.
# ==> Compare our SCF to Psi4 <==
# Call psi4.energy() to compute the SCF energy
SCF_E_psi = psi4.energy('SCF')
psi4.core.clean()
print(f"PSI4 {SCF_E_psi:.6f} Ours {SCF_E:.6f} Absolute difference {np.abs(SCF_E-SCF_E_psi)}")
Bonus Exercise 27
Modify the value of E_conv and describe its effect the number of iterations.