<jupyter_start><jupyter_text>Unrestricted Hartree-Fock**This assignment is due on Friday, January 31 at 5PM by email.**In carrying out the implementation of the closed-shell Hartree-Fock method we have, in fact, had to solve most of the design problems associated with open-shell Hartree-Fock.\begin{align}E_{elec}^{UHF} =&\ \sum_i^{n_\alpha} \ h_{ii}^{\alpha} + \sum_i^{n_\beta} \ h_{ii}^{\beta} \\               &+ \sum_{i<j}^{n_\alpha} \biggl\{ (i^\alpha i^\alpha | j^\alpha j^\alpha) - (i^\alpha j^\alpha | i^\alpha j^\alpha) \biggr\} + \sum_{i<j}^{n_\beta} \biggl\{ (i^\beta i^\beta | j^\beta j^\beta) - (i^\beta j^\beta | i^\beta j^\beta) \biggr\} \\               &+ \sum_i^{n_\alpha} \sum_j^{n_\beta} (i^\alpha i^\alpha | j^\beta j^\beta)\\               =&\ \sum_i^{n_\alpha} \ h_{ii}^{\alpha} + \sum_i^{n_\beta} \ h_{ii}^{\beta} \\               & + \sum_{i<j}^{n_\alpha} \biggl\{ J_{ij}^{\alpha} - K_{ij}^{\alpha} \biggr\} + \sum_{i<j}^{n_\beta} \biggl\{ J_{ij}^{\beta} - K_{ij}^{\beta} \biggr\} \\               & + \sum_i^{n_\alpha} \sum_j^{n_\beta} J_{ij}^{\alpha\beta} \\\end{align} The density matrix\begin{align}D_{\mu\nu}^\alpha =&\ \sum_i^{n_{\alpha}} C_{\mu i}^\alpha C_{\nu i}^\alpha \\D_{\mu\nu}^\beta =&\ \sum_i^{n_{\beta}} C_{\mu i}^\beta C_{\nu i}^\beta\end{align}  One-electron energy\begin{align}E_{one} =&\ \sum_i^{n_\alpha} h_{ii}^\alpha + \sum_i^{n_\beta} h_{ii}^\beta \\        =&\ \sum_i^{n_\alpha} \biggl\{ \sum_{\mu \nu}^{AO} C_{\mu i}^\alpha C_{\nu i}^\alpha h_{\mu \nu} \biggr\} + \sum_i^{n_\beta} \biggl\{ \sum_{\mu \nu}^{AO} C_{\mu i}^\beta C_{\nu i}^\beta h_{\mu \nu} \biggr\} \\        =&\ \sum_{\mu \nu}^{AO} \biggl\{ \sum_i^{n_\alpha} C_{\mu i}^\alpha C_{\nu i}^\alpha + C_{\mu i}^\beta C_{\nu i}^\beta \biggr\} h_{\mu \nu} \\        =&\ \sum_{\mu \nu}^{AO} \biggl\{ D_{\mu \nu}^\alpha + D_{\mu \nu}^\beta \biggr\} h_{\mu \nu}\end{align}Notice that the one-electron energy is composed of an $\alpha$ and $\beta$ energy.  Two-electron energy\begin{align}E_{two} =&\ \sum_{i<j}^{n_\alpha} \biggl\{ J_{ij}^{\alpha} - K_{ij}^{\alpha} \biggr\} + \sum_{i<j}^{n_\beta} \biggl\{ J_{ij}^{\beta} - K_{ij}^{\beta} \biggr\} + \sum_i^{n_\alpha} \sum_j^{n_\beta} J_{ij}^{\alpha\beta} + \sum_i^{n_\beta} \sum_j^{n_\alpha} J_{ij}^{\beta\alpha}  \\\\        =&\ \sum_{i<j}^{n_\alpha} \Biggl[ \sum_{\mu\nu\rho\sigma}^{AO} \biggl\{ C_{\mu i}^\alpha C_{\nu i}^\alpha C_{\rho j}^\alpha C_{\sigma j}^\alpha (\mu\nu | \rho\sigma) - C_{\mu i}^\alpha C_{\rho j}^\alpha C_{\nu i}^\alpha C_{\sigma j}^\alpha (\mu\rho | \nu\sigma)\biggr\} \Biggr] \\        +&\ \sum_{i<j}^{n_\beta} \Biggl[ \sum_{\mu\nu\rho\sigma}^{AO} \biggl\{ C_{\mu i}^\beta C_{\nu i}^\beta C_{\rho j}^\beta C_{\sigma j}^\beta (\mu\nu | \rho\sigma) - C_{\mu i}^\beta C_{\rho j}^\beta C_{\nu i}^\beta C_{\sigma j}^\beta (\mu\rho | \nu\sigma)\biggr\} \Biggr] \\        +&\ \sum_i^{n_\alpha} \sum_j^{n_\beta} \Biggl[ \sum_{\mu\nu\rho\sigma}^{AO} \biggl\{ C_{\mu i}^\alpha C_{\nu i}^\alpha C_{\rho j}^\beta C_{\sigma j}^\beta (\mu\nu | \rho\sigma) \biggr\} \Biggr] \\        +&\ \sum_i^{n_\beta} \sum_j^{n_\alpha} \Biggl[ \sum_{\mu\nu\rho\sigma}^{AO} \biggl\{ C_{\mu i}^\beta C_{\nu i}^\beta C_{\rho j}^\alpha C_{\sigma j}^\alpha (\mu\nu | \rho\sigma) \biggr\} \Biggr] \\\\        =&\ \sum_{\mu\nu\rho\sigma}^{AO} \Biggl[ \biggl\{ \sum_i^{n_\alpha} C_{\mu i}^\alpha C_{\nu i}^\alpha \biggr\} \biggl\{ \sum_j^{n_\alpha} C_{\rho j}^\alpha C_{\sigma j}^\alpha \biggr\} \biggl\{ (\mu\nu|\rho\sigma) - (\mu\rho|\nu\sigma) \biggr\} \Biggr] \\        +&\ \sum_{\mu\nu\rho\sigma}^{AO} \Biggl[ \biggl\{ \sum_i^{n_\beta} C_{\mu i}^\beta C_{\nu i}^\beta \biggr\} \biggl\{ \sum_j^{n_\beta} C_{\rho j}^\beta C_{\sigma j}^\beta \biggr\} \biggl\{ (\mu\nu|\rho\sigma) - (\mu\rho|\nu\sigma) \biggr\} \Biggr] \\        +&\ \sum_{\mu\nu\rho\sigma}^{AO} \Biggl[ \biggl\{ \sum_i^{n_\alpha} C_{\mu i}^\alpha C_{\nu i}^\alpha \biggr\} \biggl\{ \sum_j^{n_\beta} C_{\rho j}^\beta C_{\sigma j}^\beta \biggr\} (\mu\nu|\rho\sigma) \Biggr] \\        +&\ \sum_{\mu\nu\rho\sigma}^{AO} \Biggl[ \biggl\{ \sum_i^{n_\beta} C_{\mu i}^\beta C_{\nu i}^\beta \biggr\} \biggl\{ \sum_j^{n_\alpha} C_{\rho j}^\alpha C_{\sigma j}^\alpha \biggr\} (\mu\nu|\rho\sigma) \Biggr] \\\\        =&\ \sum_{\mu\nu\rho\sigma}^{AO} D_{\mu\nu}^\alpha D_{\rho\sigma}^\alpha \biggl\{ (\mu\nu|\rho\sigma) - (\mu\rho|\nu\sigma) \biggr\} \\        +&\ \sum_{\mu\nu\rho\sigma}^{AO} D_{\mu\nu}^\beta D_{\rho\sigma}^\beta \biggl\{ (\mu\nu|\rho\sigma) - (\mu\rho|\nu\sigma) \biggr\} \\        +&\ \sum_{\mu\nu\rho\sigma}^{AO} D_{\mu\nu}^\alpha D_{\rho\sigma}^\beta (\mu\nu|\rho\sigma) \\        +&\ \sum_{\mu\nu\rho\sigma}^{AO} D_{\mu\nu}^\beta D_{\rho\sigma}^\alpha (\mu\nu|\rho\sigma) \\\end{align}  Fock Matrix\begin{align}F_{\mu\nu}^{\alpha} =&\ h_{\mu\nu}^\alpha + \sum_{\rho\sigma}^{AO} D_{\rho\sigma}^\alpha \biggl\{ (\mu\nu|\rho\sigma) - (\mu\rho|\nu\sigma) \biggr\}  + D_{\rho\sigma}^\beta (\mu\nu|\rho\sigma) \\F_{\mu\nu}^{\alpha} =&\ h_{\mu\nu}^\alpha + J_{\mu\nu}^{\alpha} - K_{\mu\nu}^\alpha + J_{\mu\nu}^\beta\\\\F_{\mu\nu}^{\beta} =&\ h_{\mu\nu}^\beta + \sum_{\rho\sigma}^{AO} D_{\rho\sigma}^\beta \biggl\{ (\mu\nu|\rho\sigma) - (\mu\rho|\nu\sigma) \biggr\}  + D_{\rho\sigma}^\alpha (\mu\nu|\rho\sigma) \\F_{\mu\nu}^{\beta} =&\ h_{\mu\nu}^\beta + J_{\mu\nu}^{\beta} - K_{\mu\nu}^\beta + J_{\mu\nu}^\alpha\end{align} Electronic Energy\begin{align}E_{elec} = \frac{1}{2}\Biggl[ \sum_{\mu \nu}^{AO} \biggl\{ D_{\mu\nu}^\alpha + D_{\mu\nu}^\beta \biggr\} h_{\mu\nu} + \sum_{\mu \nu}^{AO} D_{\mu\nu}^\alpha F_{\mu\nu}^\alpha + \sum_{\mu\nu}^{AO} D_{\mu\nu}^\beta F_{\mu\nu}^{\beta} \Biggr]\end{align}The $\tfrac{1}{2}$ is to correct for double counting.Additional insights into UHF can be found in Modern Quantum Chemistry by Szabo and Ostlund at Chapter 3 Section 8. Spin Contamination Metric\begin{align}S_{i j}^{\mathrm{metric}} =& C_{\mu i}^{\alpha} S_{\mu \nu} C_{\nu j}^\beta \\S^{\mathrm{contamination}} =& \mathrm{min}(n_{\alpha}, n_{\beta}) - (S^{\mathrm{metric}} \cdot S^{\mathrm{metric}})\end{align}  ProcedureFor this assignment, please create a separate `.py` file as this will form the base for the remaining assignments.1. Load Psi4 and NumPy2. Setup your initial conditions   * Construct a new molecule in C1 symmetry   * Load in the basis set   * Create a `MintsHelper` object to compute the integrals.3. Compute needed integrals4. Form the one-electron Hamiltonian.5. Construct the orthogonalizer $S^{-1/2}$6. Construct the initial density matrices $D^\alpha$ and $D^\beta$   * Form the initial (transformed) Fock matrices using a core Hamiltonian guess   * Diagonalize the Fock matrices   * Form the SCF eigenvectors $C^\alpha$ and $C^\beta$   * Form the initial density matrices $D^\alpha$ and $D^\beta$7. Perform SCF iterations until convergence is achieved   * Construct new Fock matrices   * Calculate this iteration's electronic and total energies   * Transform the Fock matrices to the orthonormal basis   * Diagonalize the Fock matrices   * Construct the new SCF eigenvector matrices   * Form the new density matrices   * Test convergency of the density matrices and energy8. Compute and print the spin-contamination. My output:    UHF iteration   1: energy   -73.50799290557293  dE -7.35080E+01    UHF iteration   2: energy   -74.66361368523238  dE -1.15562E+00    UHF iteration   3: energy   -74.66614383044545  dE -2.53015E-03    UHF iteration   4: energy   -74.66638263555237  dE -2.38805E-04    UHF iteration   5: energy   -74.66644299107340  dE -6.03555E-05    UHF iteration   6: energy   -74.66646469386085  dE -2.17028E-05    UHF iteration   7: energy   -74.66647352191376  dE -8.82805E-06    UHF iteration   8: energy   -74.66647726970710  dE -3.74779E-06    UHF iteration   9: energy   -74.66647888623531  dE -1.61653E-06    UHF iteration  10: energy   -74.66647958779794  dE -7.01563E-07    UHF iteration  11: energy   -74.66647989301480  dE -3.05217E-07    UHF iteration  12: energy   -74.66648002592984  dE -1.32915E-07    UHF iteration  13: energy   -74.66648008383399  dE -5.79041E-08    UHF iteration  14: energy   -74.66648010906377  dE -2.52298E-08    UHF iteration  15: energy   -74.66648012005750  dE -1.09937E-08    UHF iteration  16: energy   -74.66648012484804  dE -4.79054E-09    UHF iteration  17: energy   -74.66648012693554  dE -2.08750E-09    UHF iteration  18: energy   -74.66648012784520  dE -9.09665E-10    UHF iteration  19: energy   -74.66648012824157  dE -3.96369E-10    UHF iteration  20: energy   -74.66648012841434  dE -1.72761E-10    UHF iteration  21: energy   -74.66648012848960  dE -7.52607E-11        Spin Contamination Metric: 6.40288E-03     Testing against Psi4<jupyter_code>import psi4

psi4.set_options({'basis': 'sto-3g',
                 'scf_type': 'pk',
                 'e_convergence': 1e-10,
                 'reference': 'uhf'})
psi4.energy('scf')
import numpy as np
from openshell import Settings

np.set_printoptions(suppress=True, linewidth=200)
molecule = psi4.geometry(Settings['molecule'])
molecule.update_geometry()
nalpha = Settings['nalpha']
nbeta = Settings['nbeta']

scf_max_iter = Settings['scf_max_iter']

basis = psi4.core.BasisSet.build(molecule, 'BASIS', Settings['basis'])
mints = psi4.core.MintsHelper(basis)
S = mints.ao_overlap().np
T = mints.ao_kinetic().np
V = mints.ao_potential().np
I = mints.ao_eri().np

H = T + V
print(H)
from scipy.linalg import fractional_matrix_power

A = fractional_matrix_power(S, -0.5)
print(A)
Fa = A.dot(H).dot(A)
Fb = A.dot(H).dot(A)

_, Ca = np.linalg.eigh(Fa)
_, Cb = np.linalg.eigh(Fb)

Ca = A.dot(Ca)
Cb = A.dot(Cb)

Cocca = Ca[:, :nalpha]
Coccb = Cb[:, :nbeta]

Da = np.einsum('pi,qi->pq', Cocca, Cocca)
Db = np.einsum('pi,qi->pq', Coccb, Coccb)

print(Da)
print(Db)<jupyter_output>[[ 1.06661727 -0.29755998 -0.02536219 -0.         -0.          0.04007352  0.04007352]
 [-0.29755998  1.32589965  0.15348429  0.          0.         -0.17684186 -0.17684186]
 [-0.02536219  0.15348429  1.08573619  0.          0.         -0.11682694 -0.11682694]
 [-0.          0.          0.          1.         -0.         -0.          0.        ]
 [-0.          0.          0.         -0.          1.17091368  0.1776736  -0.1776736 ]
 [ 0.04007352 -0.17684186 -0.11682694 -0.          0.1776736   0.05924138  0.00532125]
 [ 0.04007352 -0.17684186 -0.11682694  0.         -0.1776736   0.00532125  0.05924138]]
[[ 1.06661727 -0.29755998 -0.02536219 -0.         -0.          0.04007352  0.04007352]
 [-0.29755998  1.32589965  0.15348429  0.          0.         -0.17684186 -0.17684186]
 [-0.02536219  0.15348429  1.08573619 -0.          0.         -0.11682694 -0.11682694]
 [-0.          0.         -0.          0.         -0.         -0.          0.        ]
 [-0.          0.          0.         -0. [...]
{'fraction_non_alphanumeric': 0.21825286073123082, 'fraction_numerical': 0.09479951623406828, 'mean_word_length': 2.645303492709393, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '462331', 'license': 'Apache-2.0', 'n_tokens_mistral': 4805, 'n_tokens_neox': 4125, 'n_words': 930, 'path': '2020/programming/2_UHF.ipynb', 'repo_name': 'GQCG-edu/chem-8950', 'size': '10749'}