<jupyter_start><jupyter_text>This lab activity is intended to teach students the mechanics of the Hartree-Fock procedure, without getting into the details of calculating the 1 and 2 electron integrals, by using the machinary of the Psi4 quantum chemistry software package.Prerequisite knowledge:  This lab assumes that the students understand the principles, though not the mathematics, of constructing the Fock matrix and that diagaonlizing the Fock matrix leads to the molecular orbital energies and coefficients.  The lab also assumes all the standard python pre-requisites of all Psi4Education labs.  Learning ojbectives:1. Recognize that the AO basis is not orthonormal and must be transformed.2. Transform a from one basis to another using a transformation matrix.3. Recognize the interative nature of the HF procedure and what leds to the need for the iterative process.  4. Define convergence and use convergence criteria in a self-consistent procedure. 5. Calculate MO energies and coefficients by diagonalizing the Fock matrix.  Authors: <NAME> (<EMAIL>; ORCID: 0000-0002-4381-1239) and <NAME> (<EMAIL>; ORCID: )Copyright: Psi4Education Project, 2019  The Hartree-Fock procedureIn this lab activity, you will be building and diagonalizing 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.<jupyter_code># ==> Import Psi4 & NumPy <==
import psi4
import numpy as np<jupyter_output><empty_output><jupyter_text>Specifying the moleculeBefore we can begin to implement the HF procedure, we need to specifcy 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.<jupyter_code># ==> Set Basic Psi4 Options <==
# Memory specification
psi4.set_memory(int(5e8))
numpy_memory = 2

# Set output file
psi4.core.set_output_file('output.dat', False)

# Define Physicist's water -- don't forget C1 symmetry!
mol = psi4.geometry("""
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

# Set computation options
psi4.set_options({'basis': 'cc-pvdz',
                  'scf_type': 'pk',
                  'e_convergence': 1e-8})<jupyter_output><empty_output><jupyter_text>The Hartree Fock procedure is iterative (more on this later).  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.<jupyter_code># ==> Set default program options <==
# Maximum SCF iterations
MAXITER = 40
# Energy convergence criterion
E_conv = 1.0e-6<jupyter_output><empty_output><jupyter_text>Next, we need to build the wavefunction from the basis functions.  We store the wavefunction in a variable called `wfn`.  Next we have to set up a data structure, called a class, to calculate the molecular integrals.  We will call this data structure `mints`.<jupyter_code># ==> Compute static 1e- and 2e- quantities with Psi4 <==
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('basis'))
mints = psi4.core.MintsHelper(wfn.basisset())<jupyter_output><empty_output><jupyter_text>We use the function `ao_overlap` to calculate the overlap integrals between all the AO basis functions.  We cast this to a numpy array called `S`.<jupyter_code># Overlap matrix
S = np.asarray(mints.ao_overlap())<jupyter_output><empty_output><jupyter_text>Now that we have constructed a matrix of the overlap integrals, we can determine how many basis functions we have and how many occupied orbitals there are.  Since we know how many electrons our molecule has, this can be a good check to make sure our code is working correctly so far.   In the box below, use the function `nalpha()` acting on the wavefunction object we created above to determine the number of doubly occupied orbitals.  Save this answer as a variable called `ndocc`.  This should match your expectation for the molecule you chose.Then determine the number of basis functions.  Remember that your overlap matrix is of size NxN where N is the number of basis functions.  Use the numpy documentation to look up the `shape` function and use that function to determine the size of the `S` array.  Use that information to determine how many basis functions there are and save your answer in a variable called `nbf`.<jupyter_code># This is the code we would expect the student to write
ndocc = wfn.nalpha()
size_overlap = S.shape
#print('The size of S is', size_overlap)
nbf = size_overlap[0]

print(F'Number of occupied orbitals: {ndocc}')
print(F'Number of basis functions: {nbf}')<jupyter_output>Number of occupied orbitals: 5
Number of basis functions: 24<jupyter_text>Now that our basis set checks out, we can build the electron-repulsion integral tensor, which calculates the electron repulsion between the atomic orbital wavefunctions, and the core Hamiltonian.<jupyter_code># Build ERI Tensor
I = np.asarray(mints.ao_eri())

# Build core Hamiltonian
T = np.asarray(mints.ao_kinetic())
V = np.asarray(mints.ao_potential())
H = T + V<jupyter_output><empty_output><jupyter_text>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 real basis set of Gaussian wavefuctions, this may not be the case.  We need to check our the overlap integral array and see if we have an orthonormal basis. Question:  If our basis was orthonormal, what would the overlap integral array look like?  Student Answer BoxAnswer: It would be 1's along the diagonal and 0's for all other elements. Construct an array of the same size as the overlap array (`S`) that has 1's along the diagonal and 0's everywhere else.  Use the numpy documentation to identify a function to do this.  Then compare that array to the `S` array to determine if the AO basis is orthonormal.<jupyter_code>## This is the code we would expect the student to write.
size_S = S.shape[0]
comparison_array = np.eye(size_S)
#print(comparison_array)
orthonormal_check = np.allclose(S, comparison_array)
#print(orthonormal_check)
print(F'The AO basis is orthonormal? {orthonormal_check}')<jupyter_output>The AO basis is orthonormal? False<jupyter_text>Since our AO basis set was not orthonormal, we must construct an orthogonalization matrix, `A` such that ${\bf A}^{\dagger}{\bf SA} = {\bf 1}$, and cast it to a numpy array.<jupyter_code># ==> Construct AO orthogonalization matrix A <==
A = mints.ao_overlap()
A.power(-0.5, 1.e-16)
A = np.asarray(A)<jupyter_output><empty_output><jupyter_text>Use the orthogonalization matrix `A` to transform the overlap matrix, `S`. Check the transformed overlap matrix to make sure it represents an orthonormal basis.<jupyter_code># This is the code we would expect the student to write
S_p = A.dot(S).dot(A)
transformation_check = np.allclose(S_p, comparison_array)
#print(transformation_check)

if transformation_check is True:
    print('Transformation success!')
else:
    print("Whoops...something went wrong. Check that you've correctly built the transformation matrix.")<jupyter_output>Transformation success!<jupyter_text>Now we would have to recompute the ERI and all the core Hamiltonian matrices in the new orthogonal basis **OR** we can just calculate orbital coefficient matrix in the transformed basis by diagonalizing the transformed Fock matrix, and then just go backwards and get the cofficients in the original AO basis.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.  So we start with a guess for the Fock matrix, which is the core Hamiltonian matrix that we already calculated.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.  Look up the `linalg.eigh` function in the numpy documentation, and use it to diagonalize the transformed Fock matrix.<jupyter_code># This is the code we would expect the student to write

# Transformed Fock matrix
F_p = A.dot(H).dot(A)

# Diagonalize F_p for eigenvalues & eigenvectors with NumPy
e, C_p = np.linalg.eigh(F_p)<jupyter_output><empty_output><jupyter_text>Now that we have the coefficents in the transformed basis, we need to go back and get the coefficients in the original AO basis.<jupyter_code># Transform C_p back into AO basis
C = A.dot(C_p)<jupyter_output><empty_output><jupyter_text>Next, 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.<jupyter_code># Grab occupied orbitals
C_occ = C[:, :ndocc]<jupyter_output><empty_output><jupyter_text>To build the density matrix, we will use the numpy function `einsum`, 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.  A good tutorial about `einsum` can be found at http://ajcr.net/Basic-guide-to-einsum/.  To specify the operations you want to do, you use the **Ein**stein **Sum**mation convention.<jupyter_code># Build density matrix from occupied orbitals
D = np.einsum('pi,qi->pq', C_occ, C_occ)<jupyter_output><empty_output><jupyter_text>Now we will use the density matrix to build the Fock matrix.  The code block below sets up the skelton of the Hartree-Fock procedure.  The basic steps are:1. Calculate the Fock Matrix from the density matrix and the nuclear repulsion energy (which we only need to calculate once at the beginning because we are invoking the the Born-Oppenheimer approximation thus making the nucleus stationary.  2. Calculate the energy from the Fock matrix.3. 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.  (You will write this section.)4. 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.  (You will write this section.)Follow the prompts below to complete your HF code.<jupyter_code># ==> Nuclear Repulsion Energy <==
E_nuc = mol.nuclear_repulsion_energy()

# ==> SCF Iterations <==
# Pre-iteration energy declarations
SCF_E = 0.0
E_old = 0.0

print('==> Starting SCF Iterations <==\n')

# Begin Iterations
for scf_iter in range(1, MAXITER + 1):
    # Build Fock matrix
    J = np.einsum('pqrs,rs->pq', I, D)
    K = np.einsum('prqs,rs->pq', I, D)
    F = H + 2*J - K
    
    # Compute RHF energy; this use of einsum is matrix multiplication
    SCF_E = np.einsum('pq,pq->', (H + F), D) + E_nuc
    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
    
    ## This is the code we would expect the student to write
    if (abs(SCF_E - E_old) < E_conv):
        break
    E_old = SCF_E
    
    # Compute new orbital guess.  Remember the steps:
    # 1. Transform the Fock matrix
    # 2. Diagonalize the Fock matrix.
    # 3. Tranform the coefficient matrix back to original AO basis set
    # 4. Take a slice of the coefficient matrix of just the occupied orbitals.
    # 5. Use the occupied orbitals to construct the density
    
    ## This is the code we would expect the student to write
    F_p =  A.dot(F).dot(A)
    e, C_p = np.linalg.eigh(F_p)
    C = A.dot(C_p)
    C_occ = C[:, :ndocc]
    D = np.einsum('pi,qi->pq', C_occ, C_occ)
    
    # 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]')<jupyter_output>==> Starting SCF Iterations <==

SCF Iteration 1: Energy = -68.98003273 dE = -68.98003273
SCF Iteration 2: Energy = -69.64725443 dE = -0.66722169
SCF Iteration 3: Energy = -72.84030309 dE = -3.19304867
SCF Iteration 4: Energy = -72.89488391 dE = -0.05458081
SCF Iteration 5: Energy = -74.12078065 dE = -1.22589674
SCF Iteration 6: Energy = -74.86718195 dE = -0.74640130
SCF Iteration 7: Energy = -75.41490878 dE = -0.54772683
SCF Iteration 8: Energy = -75.70767221 dE = -0.29276343
SCF Iteration 9: Energy = -75.86052421 dE = -0.15285200
SCF Iteration 10: Energy = -75.93143852 dE = -0.07091432
SCF Iteration 11: Energy = -75.96402440 dE = -0.03258587
SCF Iteration 12: Energy = -75.97842921 dE = -0.01440481
SCF Iteration 13: Energy = -75.98481699 dE = -0.00638779
SCF Iteration 14: Energy = -75.98761283 dE = -0.00279584
SCF Iteration 15: Energy = -75.98884089 dE = -0.00122806
SCF Iteration 16: Energy = -75.98937775 dE = -0.00053686
SCF Iteration 17: Energy = -75.98961294 dE = -0.00023518
SCF It[...]<jupyter_text>Use the `psi4.energy` function to have psi4 calcluate the energy of your molecule.  Compare your results.<jupyter_code># Compare to Psi4
SCF_E_psi = psi4.energy('SCF')
print(SCF_E_psi)<jupyter_output>-75.9897957853276
{'fraction_non_alphanumeric': 0.07503111046043481, 'fraction_numerical': 0.03433130810335993, 'mean_word_length': 4.484544359694902, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 78, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '72029', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4446, 'n_tokens_neox': 4015, 'n_words': 1821, 'path': 'advanced-labs/Hartree-Fock/HF_lab.ipynb', 'repo_name': 'dbmagers/labs', 'size': '13718'}