{
"cells": [
{
"cell_type": "markdown",
"id": "baada2eb",
"metadata": {},
"source": [
"# Hartree–Fock theory\n",
"\n",
"With the state of the system being described by a single Slater determinant, the Hartree–Fock wave function is given as that which minimizes the electronic energy in a variational sense with respect to variations in the spin orbitals. It represents a cornerstone in quantum chemistry and provides total electronic energies that are within 1% of the exact results and a wide range of molecular properties that are within 5–10% accuracy. Moreover, the Hartree–Fock method serves as starting points for the formulation of many other, more accurate, wave function methods as well as the Kohn–Sham formulation of density functional theory.\n",
"\n",
"(hartree-fock-equation)=\n",
"## Hartree–Fock equation\n",
"\n",
"In the Hartree–Fock approximation, the many-electron wave function takes the form of a [Slater determinant](sec:slater)\n",
"\n",
"$$\n",
" | \\Phi \\rangle =\n",
" \\frac{1}{\\sqrt{N!}}\n",
" \\begin{vmatrix}\n",
" \\psi_{1}(\\mathbf{r}_1) & \\cdots & \\psi_{N}(\\mathbf{r}_1) \\\\\n",
" \\vdots & \\ddots & \\vdots \\\\\n",
" \\psi_{1}(\\mathbf{r}_N) & \\cdots & \\psi_{N}(\\mathbf{r}_N) \\\\\n",
" \\end{vmatrix} \n",
"$$\n",
"\n",
"where $\\psi_i$ are the single-electron wave functins known as [spin orbitals](sec:orbitals). The Hartree–Fock energy and the associated state is found by minimizing the energy functional\n",
"\n",
"$$\n",
"E_\\mathrm{HF} =\n",
"\\min_{\\psi} E[\\psi]\n",
"$$\n",
"\n",
"under the constraint that the spin orbitals remain orthonormal. Here, $\\psi$ collectively refers to the entire set of $N$ spin orbitals. Such a contrained minimization is conveniently performed by means of the technique of Lagrange multipliers.\n",
"\n",
"(lagrangian)=\n",
"### Lagrangian\n",
"\n",
"In Hartree–Fock theory, we introduce the real-valued Lagrangian\n",
"\n",
"$$\n",
"L[\\psi] = E[\\psi] - \\sum_{i,j=1}^N\n",
"\\varepsilon_{ji} \\big(\n",
"\\langle \\psi_i | \\psi_j \\rangle - \\delta_{ij}\n",
"\\big)\n",
"$$\n",
"\n",
"and search for the set of spin orbitals, $\\psi$, that results in a first variation that vanishes\n",
"\n",
"$$\n",
"\\delta L = 0\n",
"$$\n",
"\n",
"Expressing the energy as the expectation value of the electronic Hamiltonian with respect to a Slater determinant and using the general expressions for [matrix elements](sec:matrix-elements), we arrive at\n",
"\n",
"\\begin{align*}\n",
"\\delta L & =\n",
"\\sum_{i=1}^N\n",
"\\langle \\delta \\psi_i | \\hat{h} | \\psi_i \\rangle +\n",
"\\sum_{i,j=1}^N\n",
"\\big(\n",
"\\langle \\delta \\psi_i \\psi_j | \\hat{g} | \\psi_i \\psi_j\\rangle -\n",
"\\langle \\delta \\psi_i \\psi_j | \\hat{g} | \\psi_j \\psi_i\\rangle\n",
"-\n",
"\\varepsilon_{ji}\n",
"\\langle \\delta \\psi_i | \\psi_j \\rangle\n",
"\\big) +\n",
"\\mbox{complex conjugate} \\\\\n",
"&=\n",
"\\sum_{i=1}^N\n",
"\\langle \\delta \\psi_i | \\big( \n",
"\\hat{f} | \\psi_i \\rangle -\n",
"\\sum_{j=1}^N\n",
"\\varepsilon_{ji} | \\psi_j \\rangle\n",
"\\big) +\n",
"\\mbox{complex conjugate}\n",
"\\end{align*}\n",
"\n",
"where we have introduced the *one-electron* Fock operator\n",
"\n",
"$$\n",
"\\hat{f} = \\hat{h} + \\sum_{j=1}^N \\big( \\hat{J}_j - \\hat{K}_j \\big)\n",
"$$\n",
"\n",
"with\n",
"\n",
"\\begin{align*}\n",
"\\hat{J}_j | \\psi_i \\rangle & = \n",
"\\Big[ \n",
"\\int \n",
"\\frac{e^2 |\\psi_j(\\mathbf{r}')|^2}{4\\pi\\varepsilon_0 |\\mathbf{r} - \\mathbf{r}'|}\n",
"d^3\\mathbf{r}'\n",
"\\Big]\n",
"| \\psi_i \\rangle \\\\\n",
"%\n",
"\\hat{K}_j | \\psi_i \\rangle & = \n",
"\\Big[\n",
"\\int\n",
"\\frac{e^2 \\psi_j^\\dagger(\\mathbf{r}')\\psi_i(\\mathbf{r}')}{4\\pi\\varepsilon_0 |\\mathbf{r} - \\mathbf{r}'|}\n",
"d^3\\mathbf{r}'\n",
"\\Big]\n",
"| \\psi_j \\rangle \n",
"\\end{align*}\n",
"\n",
"Since the first-order variation in the Lagrangian is required to vanish for general variations in the spin orbitals, we have shown that the Hartree–Fock solution is given by \n",
"\n",
"$$\n",
"\\hat{f} | \\psi_i \\rangle -\n",
"\\sum_{j=1}^N\n",
"\\varepsilon_{ji} | \\psi_j \\rangle = 0\n",
"$$\n",
"\n",
"This equation is known as the Hartree–Fock equation and it to be solved for the spin orbitals and the associated Lagrange multipliers. We note that the matrix elements of the Fock operator equal the multipliers\n",
"\n",
"$$\n",
"f_{ki} =\n",
"\\langle \\psi_k | \\hat{f} | \\psi_i \\rangle = \n",
"\\sum_{j=1}^N\n",
"\\varepsilon_{ji} \\langle \\psi_k | \\psi_j \\rangle = \\varepsilon_{ki}\n",
"$$\n",
"\n",
"(sec:canonical_hf)=\n",
"### Canonical form\n",
"Apart from a trivial overall phase factor, [unitary transformations among the occupied orbitals](sec:unitary) are shown to leave the Hartree–Fock wave function unchanged. We introduce a unitary transformation that diagonalizes the Hermitian Fock matrix\n",
"\n",
"$$\n",
"f' = \\langle \\overline{\\psi}' | \\hat{f} | \\overline{\\psi}' \\rangle \n",
"= U^\\dagger \\langle \\overline{\\psi} | \\hat{f} | \\overline{\\psi} \\rangle U =\n",
"U^\\dagger f U \n",
"$$\n",
"\n",
"We have here adopted the compact [overline notation](sec:lcao) of orbitals. In this basis of *canonical spin orbitals*, the Hartree--Fock equation takes the form\n",
"\n",
"$$\n",
"\\hat{f} | \\psi_i \\rangle =\n",
"\\varepsilon_{i} | \\psi_i \\rangle\n",
"$$\n",
"\n",
"which we recognize as an eigenvalue equation introducing the *orbital energies*, $\\varepsilon_{i}$, as the eigenvalues of the Fock operator. With an infinite number of solutions to the Hartree–Fock equation, the Hartree–Fock ground state is given by employing the $N$ spin orbitals with lowest orbital energies in the Slater determinant.\n",
"\n",
"#### In AO basis\n",
"The spatial parts of the spin orbitals, or molecular orbitals (MOs), are expanded as linear combination of atomic orbitals ([LCAO](sec:lcao)). In the basis of spin atomic orbitals, the Fock matrix becomes block diagonal\n",
"\n",
"$$\n",
"\\mathbf{F} =\n",
"\\begin{pmatrix}\n",
"\\mathbf{F}^{\\alpha\\alpha} & \\mathbf{0} \\\\\n",
"\\mathbf{0} & \\mathbf{F}^{\\beta\\beta}\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"Using the [bar notation](sec:orbitals) to distinguish $\\alpha$- and $\\beta$-spin atomic orbitals, we get\n",
"\n",
"\\begin{align*}\n",
"F_{\\mu\\nu} & = F^{\\alpha\\alpha}_{\\mu\\nu} =\n",
"h_{\\mu\\nu} + \\sum_{\\gamma\\delta} \\Big(\n",
"D_{\\gamma\\delta}(\\mu\\nu|\\gamma\\delta) -\n",
"D^\\alpha_{\\gamma\\delta}(\\mu\\delta|\\gamma\\nu)\n",
"\\Big)\n",
"\\\\\n",
"F_{\\bar{\\mu}\\bar{\\nu}} & = F^{\\beta\\beta}_{\\mu\\nu} =\n",
"h_{\\mu\\nu} + \\sum_{\\gamma\\delta} \\Big(\n",
"D_{\\gamma\\delta}(\\mu\\nu|\\gamma\\delta) -\n",
"D^\\beta_{\\gamma\\delta}(\\mu\\delta|\\gamma\\nu)\n",
"\\Big)\n",
"\\\\\n",
"F_{\\mu\\bar{\\nu}} & = F_{\\bar{\\mu}\\nu} = 0\n",
"\\end{align*}\n",
"\n",
"where\n",
"\n",
"\\begin{align*}\n",
"D_{\\gamma\\delta} &= D^\\alpha_{\\gamma\\delta} + D^\\beta_{\\gamma\\delta} \\\\\n",
"D^\\alpha_{\\gamma\\delta}& =\n",
"\\sum_{j=1}^{N_\\alpha} \\big[c_{\\gamma j}^\\alpha\\big]^* c_{\\delta j}^\\alpha \n",
"; \\quad\n",
"D^\\beta_{\\gamma\\delta} =\n",
"\\sum_{j=1}^{N_\\beta} \\big[c_{\\gamma j}^\\beta\\big]^* c_{\\delta j}^\\beta \n",
"\\\\\n",
"\\end{align*}\n",
"\n",
"The canonical Hartree--Fock equation thereby takes the form\n",
"\n",
"$$\n",
"\\mathbf{F C} = \\mathbf{S C} \\boldsymbol{\\varepsilon} \\, ,\n",
"$$\n",
"\n",
"where $\\mathbf{S}$ is the [overlap matrix](sec:overlap) and $\\boldsymbol{\\varepsilon}$ is a diagonal matrix collecting the orbital energies.\n",
"\n",
"## Hartree–Fock energy\n",
"\n",
"For a given density $\\mathbf{D}$, the Hartree–Fock energy becomes equal to\n",
"\n",
"$$\n",
"E_\\mathrm{HF} =\n",
"\\frac{1}{2}\n",
"\\mathrm{tr}\n",
"\\big[ (\\mathbf{h} + \\mathbf{F}) \\mathbf{D} \\big] + V^\\mathrm{n-n} \\, ,\n",
"$$\n",
"\n",
"where $V^\\mathrm{n-n}$ is the nuclear repulsion energy.\n",
"\n",
"## Koopmans' theorem\n",
"\n",
"The orbital energies of occupied and unoccupied orbitals, respectively, equal\n",
"\n",
"\\begin{align*}\n",
"\\varepsilon_i & = \\langle \\psi_i |\\hat{f} | \\psi_i \\rangle =\n",
"\\langle \\psi_i |\\hat{h} | \\psi_i \\rangle +\n",
"\\sum_{j\\neq i}^N\n",
"\\big(\n",
"\\langle \\psi_i | \\hat{J}_j | \\psi_i \\rangle -\n",
"\\langle \\psi_i | \\hat{K}_j | \\psi_i \\rangle \n",
"\\big) \\\\\n",
"\\varepsilon_a & = \\langle \\psi_a |\\hat{f} | \\psi_a \\rangle =\n",
"\\langle \\psi_a |\\hat{h} | \\psi_a \\rangle +\n",
"\\sum_{j=1}^N\n",
"\\big(\n",
"\\langle \\psi_a | \\hat{J}_j | \\psi_a \\rangle -\n",
"\\langle \\psi_a | \\hat{K}_j | \\psi_a \\rangle \n",
"\\big)\n",
"\\end{align*}\n",
"\n",
"where the cancellation between Coulomb and exchange terms for $j=i$ has been used in the former case. It thus appears as if $\\varepsilon_i$ relates to the energy of an electron interacting with $(N-1)$ other electrons, whereas $\\varepsilon_a$ relates to the energy of an electron interacting with $N$ other electrons. In accordance with these observations, it is readily shown from the expressions for [matrix elements](sec:matrix-elements) that\n",
"\n",
"\\begin{align*}\n",
"\\mathrm{IP} &= \n",
"E_i^{N-1} - E_\\mathrm{HF}^N = - \\varepsilon_i \\qquad (\\textsf{ionization potential}) \\\\\n",
"\\mathrm{EA} &= \n",
"E_\\mathrm{HF}^N - E_a^{N+1} = - \\varepsilon_a \\qquad (\\textsf{electron affinity})\n",
"\\end{align*}\n",
"\n",
"where, in the frozen orbital approximation, $E_i^{N-1}$ is the energy of the system after the removal of the electron in spin orbital $i$ and $E_a^{N+1}$ is the energy of the system after the addition of an electron in spin orbital $s$.\n",
"\n",
"## Brillouin's theorem\n",
"\n",
"Based on the expressions for [matrix elements](sec:matrix-elements), we find\n",
"\n",
"$$\n",
"\\langle \\Phi_\\mathrm{HF} | \\hat{H} | \\Phi_i^a \\rangle = \n",
"\\langle \\psi_i | \\hat{f} | \\psi_a \\rangle = 0 \n",
"$$\n",
"\n",
"which shows that there is no coupling between the Hartree--Fock ground state and singly excited determinants. This result is known as the *Brillouin theorem*. \n",
"\n",
"(scf-procedure)=\n",
"## SCF procedure\n",
"\n",
"Due to the summation over occupied spin orbitals that expresses the effective electron interactions, the Fock operator depends on its eigenfunctions and the canonical Hartree--Fock equation is therefore solved iteratively by means of a self-consistent field (SCF) procedure as illustrated in {numref}`rh-scf-fig`.\n",
"\n",
"```{figure} ../../img/misc/rh-scf.*\n",
"---\n",
"name: rh-scf-fig\n",
"---\n",
"Self-consistent field solution of the Hatree–Fock equation by means of the Rothaan–Hall algorithm.\n",
"```\n",
"\n",
"### Rothaan–Hall scheme\n",
"\n",
"In the following, we will consider the spin-restricted formulation where $\\alpha$- and $\\beta$-spin orbitals have identical spatial parts. We also restrict the situation to the common case of a closed-shell system such that\n",
"\n",
"\\begin{align*}\n",
"N_\\alpha & = N_\\beta = \\frac{1}{2} N \\\\\n",
"D^\\alpha_{\\gamma\\delta} & = D^\\beta_{\\gamma\\delta} = \\frac{1}{2} D_{\\gamma\\delta} =\n",
"\\sum_{j=1}^{N/2} c_{\\gamma j}^* c_{\\delta j}\n",
"\\end{align*}\n",
"\n",
"```{note}\n",
"When referring to closed-shell systems it is customary to refer to the density matrix as that for either of the spin components. \n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "dc6223d1",
"metadata": {
"tags": [
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"* Warning * Environment variable OMP_NUM_THREADS not set.\n",
"* Warning * Setting OMP_NUM_THREADS to 4.\n"
]
}
],
"source": [
"import numpy as np\n",
"import veloxchem as vlx"
]
},
{
"cell_type": "markdown",
"id": "8ec504d8",
"metadata": {},
"source": [
"# Writing your own SCF program\n",
"\n",
"Let us implement the presented Hartree–Fock SCF procedure. We will use the water molecule as an example and, as reference, we will first compute the Hartree–Fock energy using the built-in `compute` method in VeloxChem. \n",
"\n",
"## Setting up the system"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ed8acd4e-cb75-4f33-8978-ab05c446c386",
"metadata": {
"tags": [
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"* Info * Reading basis set from file: /Users/panor/opt/miniconda3/envs/vlxnb/lib/python3.10/site-packages/veloxchem/basis/CC-PVDZ\n",
" \n",
" Molecular Basis (Atomic Basis) \n",
" ================================ \n",
" \n",
" Basis: CC-PVDZ \n",
" \n",
" Atom Contracted GTOs Primitive GTOs \n",
" \n",
" O (3S,2P,1D) (17S,4P,1D) \n",
" H (2S,1P) (4S,1P) \n",
" \n",
" Contracted Basis Functions : 24 \n",
" Primitive Basis Functions : 48 \n",
" \n"
]
}
],
"source": [
"mol_str = \"\"\"\n",
"O 0.000000000000 0.000000000000 0.000000000000\n",
"H 0.000000000000 0.740848095288 0.582094932012\n",
"H 0.000000000000 -0.740848095288 0.582094932012\n",
"\"\"\"\n",
"\n",
"molecule = vlx.Molecule.read_str(mol_str, units=\"angstrom\")\n",
"basis = vlx.MolecularBasis.read(molecule, \"cc-pVDZ\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "17b35158-0505-4223-a09a-32fd6afb073d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of contracted basis functions: 24\n",
"Number of doubly occupied molecular orbitals: 5\n",
"Nuclear repulsion energy (in a.u.): 9.343638157670\n"
]
}
],
"source": [
"norb = basis.get_dimensions_of_basis(molecule)\n",
"nocc = molecule.number_of_electrons() // 2\n",
"V_nuc = molecule.nuclear_repulsion_energy()\n",
"\n",
"print(\"Number of contracted basis functions:\", norb)\n",
"print(\"Number of doubly occupied molecular orbitals:\", nocc)\n",
"print(f\"Nuclear repulsion energy (in a.u.): {V_nuc : 14.12f}\")"
]
},
{
"cell_type": "markdown",
"id": "15ef811d",
"metadata": {},
"source": [
"## Reference calculation\n",
"\n",
"Let us first perform an reference calculation using the restricted closed-shell SCF driver in VeloxChem."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "43273cce",
"metadata": {
"scrolled": true,
"tags": [
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
" Self Consistent Field Driver Setup \n",
" ==================================== \n",
" \n",
" Wave Function Model : Spin-Restricted Hartree-Fock \n",
" Initial Guess Model : Superposition of Atomic Densities \n",
" Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace \n",
" Max. Number of Iterations : 50 \n",
" Max. Number of Error Vectors : 10 \n",
" Convergence Threshold : 1.0e-06 \n",
" ERI Screening Scheme : Cauchy Schwarz + Density \n",
" ERI Screening Mode : Dynamic \n",
" ERI Screening Threshold : 1.0e-12 \n",
" Linear Dependence Threshold : 1.0e-06 \n",
" \n",
"* Info * Nuclear repulsion energy: 9.3436381577 a.u. \n",
" \n",
"* Info * Overlap matrix computed in 0.00 sec. \n",
" \n",
"* Info * Kinetic energy matrix computed in 0.00 sec. \n",
" \n",
"* Info * Nuclear potential matrix computed in 0.00 sec. \n",
" \n",
"* Info * Orthogonalization matrix computed in 0.00 sec. \n",
" \n",
"* Info * SAD initial guess computed in 0.00 sec. \n",
" \n",
"* Info * Starting Reduced Basis SCF calculation... \n",
"* Info * ...done. SCF energy in reduced basis set: -75.979046359571 a.u. Time: 0.09 sec. \n",
" \n",
"* Info * Overlap matrix computed in 0.00 sec. \n",
" \n",
"* Info * Kinetic energy matrix computed in 0.00 sec. \n",
" \n",
"* Info * Nuclear potential matrix computed in 0.00 sec. \n",
" \n",
"* Info * Orthogonalization matrix computed in 0.00 sec. \n",
" \n",
" \n",
" Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change \n",
" -------------------------------------------------------------------------------------------- \n",
" 1 -76.025869744700 0.0000000000 0.09892066 0.01121867 0.00000000 \n",
" 2 -76.026932827150 -0.0010630825 0.01920269 0.00294570 0.02594960 \n",
" 3 -76.026981056596 -0.0000482294 0.00412658 0.00076652 0.00507499 \n",
" 4 -76.026983601529 -0.0000025449 0.00247905 0.00043111 0.00166621 \n",
" 5 -76.026984175646 -0.0000005741 0.00021766 0.00003639 0.00052311 \n",
" 6 -76.026984187024 -0.0000000114 0.00003206 0.00000540 0.00012084 \n",
" 7 -76.026984187255 -0.0000000002 0.00000208 0.00000029 0.00001663 \n",
" 8 -76.026984187256 -0.0000000000 0.00000053 0.00000008 0.00000087 \n",
" \n",
" *** SCF converged in 8 iterations. Time: 0.94 sec. \n",
" \n",
" Spin-Restricted Hartree-Fock: \n",
" ----------------------------- \n",
" Total Energy : -76.0269841873 a.u. \n",
" Electronic Energy : -85.3706223449 a.u. \n",
" Nuclear Repulsion Energy : 9.3436381577 a.u. \n",
" ------------------------------------ \n",
" Gradient Norm : 0.0000005316 a.u. \n",
" \n",
" \n",
" Ground State Information \n",
" ------------------------ \n",
" Charge of Molecule : 0.0 \n",
" Multiplicity (2S+1) : 1.0 \n",
" Magnetic Quantum Number (M_S) : 0.0 \n",
" \n",
" \n",
" Spin Restricted Orbitals \n",
" ------------------------ \n",
" \n",
" Molecular Orbital No. 1: \n",
" -------------------------- \n",
" Occupation: 2.000 Energy: -20.54819 a.u. \n",
" ( 1 O 1s : -1.00) \n",
" \n",
" Molecular Orbital No. 2: \n",
" -------------------------- \n",
" Occupation: 2.000 Energy: -1.34520 a.u. \n",
" ( 1 O 2s : -0.44) ( 1 O 3s : -0.36) ( 2 H 1s : -0.20) \n",
" ( 3 H 1s : -0.20) \n",
" \n",
" Molecular Orbital No. 3: \n",
" -------------------------- \n",
" Occupation: 2.000 Energy: -0.70585 a.u. \n",
" ( 1 O 1p-1: 0.49) ( 1 O 2p-1: 0.22) ( 2 H 1s : 0.33) \n",
" ( 3 H 1s : -0.33) \n",
" \n",
" Molecular Orbital No. 4: \n",
" -------------------------- \n",
" Occupation: 2.000 Energy: -0.57109 a.u. \n",
" ( 1 O 2s : -0.15) ( 1 O 3s : -0.36) ( 1 O 1p0 : 0.55) \n",
" ( 1 O 2p0 : 0.36) ( 2 H 1s : 0.21) ( 3 H 1s : 0.21) \n",
" \n",
" Molecular Orbital No. 5: \n",
" -------------------------- \n",
" Occupation: 2.000 Energy: -0.49457 a.u. \n",
" ( 1 O 1p+1: -0.63) ( 1 O 2p+1: -0.49) \n",
" \n",
" Molecular Orbital No. 6: \n",
" -------------------------- \n",
" Occupation: 0.000 Energy: 0.18787 a.u. \n",
" ( 1 O 3s : 1.02) ( 1 O 1p0 : 0.19) ( 1 O 2p0 : 0.33) \n",
" ( 2 H 2s : -0.83) ( 3 H 2s : -0.83) \n",
" \n",
" Molecular Orbital No. 7: \n",
" -------------------------- \n",
" Occupation: 0.000 Energy: 0.25852 a.u. \n",
" ( 1 O 1p-1: 0.28) ( 1 O 2p-1: 0.67) ( 2 H 2s : -1.48) \n",
" ( 3 H 2s : 1.48) \n",
" \n",
" Molecular Orbital No. 8: \n",
" -------------------------- \n",
" Occupation: 0.000 Energy: 0.79749 a.u. \n",
" ( 1 O 1p-1: -0.26) ( 1 O 2p-1: -0.49) ( 2 H 1s : 0.95) \n",
" ( 2 H 2s : -0.66) ( 2 H 1p0 : 0.16) ( 3 H 1s : -0.95) \n",
" ( 3 H 2s : 0.66) ( 3 H 1p0 : -0.16) \n",
" \n",
" Molecular Orbital No. 9: \n",
" -------------------------- \n",
" Occupation: 0.000 Energy: 0.87271 a.u. \n",
" ( 1 O 2s : 0.25) ( 1 O 3s : -0.34) ( 1 O 1p0 : 0.34) \n",
" ( 2 H 1s : -0.77) ( 2 H 2s : 0.54) ( 2 H 1p-1: -0.31) \n",
" ( 3 H 1s : -0.77) ( 3 H 2s : 0.54) ( 3 H 1p-1: 0.31) \n",
" \n",
" Molecular Orbital No. 10: \n",
" -------------------------- \n",
" Occupation: 0.000 Energy: 1.16315 a.u. \n",
" ( 1 O 3s : 0.78) ( 1 O 1p0 : -0.74) ( 1 O 2p0 : 1.31) \n",
" ( 2 H 1s : -0.59) ( 2 H 1p0 : 0.24) ( 3 H 1s : -0.59) \n",
" ( 3 H 1p0 : 0.24) \n",
" \n"
]
}
],
"source": [
"scf_drv = vlx.ScfRestrictedDriver()\n",
"scf_drv.compute(molecule, basis)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "b7a9c77c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Final HF energy: -76.02698419 Hartree\n"
]
}
],
"source": [
"print(f\"Final HF energy: {scf_drv.get_scf_energy() : 12.8f} Hartree\")"
]
},
{
"cell_type": "markdown",
"id": "c8165a35",
"metadata": {},
"source": [
"## Getting integrals in AO basis"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "c411f050",
"metadata": {},
"outputs": [],
"source": [
"# overlap matrix\n",
"overlap_drv = vlx.OverlapIntegralsDriver()\n",
"S = overlap_drv.compute(molecule, basis).to_numpy()\n",
"\n",
"# one-electron Hamiltonian\n",
"kinetic_drv = vlx.KineticEnergyIntegralsDriver()\n",
"T = kinetic_drv.compute(molecule, basis).to_numpy()\n",
"\n",
"nucpot_drv = vlx.NuclearPotentialIntegralsDriver()\n",
"V = -nucpot_drv.compute(molecule, basis).to_numpy()\n",
"\n",
"h = T + V\n",
"\n",
"# two-electron Hamiltonian\n",
"eri_drv = vlx.ElectronRepulsionIntegralsDriver()\n",
"g = np.zeros((norb, norb, norb, norb))\n",
"eri_drv.compute_in_memory(molecule, basis, g)"
]
},
{
"cell_type": "markdown",
"id": "d27658ba",
"metadata": {},
"source": [
"## Orthogonalization of the AO basis\n",
"In order to use the `np.linalg.eigh` function from NumPy to solve the Hartree--Fock equation, we first retrieve an *orthogonal* AO (OAO) basis by means of a non-unitary transformation matrix $\\mathbf{X}$ such that\n",
"\n",
"$$\n",
"| \\overline{\\chi^\\mathrm{OAO}} \\rangle = | \\overline{\\chi} \\rangle \\mathbf{X}\n",
"$$\n",
"\n",
"and \n",
"\n",
"$$\n",
"\\mathbf{S}^\\mathrm{OAO} = \\mathbf{X}^\\dagger \\mathbf{S X} = \\mathbf{I}\n",
"$$\n",
"\n",
"In the absence of linear dependencies in the AO basis, the overlap matrix is symmetric (and positive-definite) and it can first be diagonalized by a unitary transformation\n",
"\n",
"$$\n",
"\\mathbf{U}^{\\dagger} \\mathbf{S U} = \\boldsymbol{\\sigma}\n",
"$$\n",
"\n",
"where $\\boldsymbol{\\sigma}$ is a diagonal matrix collecting the eigenvalues. It is then straightfoward to contruct explicit forms of $\\mathbf{X}$. There exist two common choices\n",
"\n",
"\\begin{align*}\n",
"\\mathbf{X} = \\mathbf{U} \\boldsymbol{\\sigma}^{-\\frac{1}{2}} \\mathbf{U}^\\dagger \\qquad & \\textsf{symmetric, or Löwdin} \\\\\n",
"\\mathbf{X} = \\mathbf{U} \\boldsymbol{\\sigma}^{-\\frac{1}{2}} \\qquad & \\textsf{canonical}\n",
"\\end{align*}\n",
"\n",
"We can readily verify that both expression satisfy $\\mathbf{X}^\\dagger \\mathbf{S X} = \\mathbf{I}$. The expression for the associated transformation of MO coefficients is determined from the relation\n",
"\n",
"$$\n",
"| \\overline{\\phi} \\rangle = | \\overline{\\chi} \\rangle \\mathbf{C} =\n",
"| \\overline{\\chi^\\mathrm{OAO}} \\rangle \\mathbf{X}^{-1} \\mathbf{C}\n",
"$$\n",
"\n",
"We identify\n",
"\n",
"$$\n",
"\\mathbf{C}^\\mathrm{OAO} = \\mathbf{X}^{-1} \\mathbf{C}\n",
"$$\n",
"\n",
"or\n",
"\n",
"$$\n",
"\\mathbf{C} = \\mathbf{X C}^\\mathrm{OAO}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "11fff6cc",
"metadata": {},
"outputs": [],
"source": [
"# symmetric transformation\n",
"sigma, U = np.linalg.eigh(S)\n",
"X = np.einsum(\"ik,k,jk->ij\", U, 1 / np.sqrt(sigma), U)"
]
},
{
"cell_type": "markdown",
"id": "96b0b7e6",
"metadata": {},
"source": [
"## Solving the Hartree–Fock equation\n",
"For a given Fock matrix, we solve the Hartree–Fock equation by the following steps:\n",
"\n",
"1. transform the Fock matrix to the OAO basis\n",
"2. diagonalize the Fock matrix\n",
"3. transform the MO coefficients back to AO basis"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "cd39c323",
"metadata": {},
"outputs": [],
"source": [
"def get_MO_coeff(F):\n",
"\n",
" F_OAO = np.einsum(\"ki,kl,lj->ij\", X, F, X)\n",
" epsilon, C_OAO = np.linalg.eigh(F_OAO)\n",
" C = np.einsum(\"ik,kj->ij\", X, C_OAO)\n",
"\n",
" return C"
]
},
{
"cell_type": "markdown",
"id": "697448c3",
"metadata": {},
"source": [
"## SCF iterations\n",
"We form an initial guess for the density based on the core Hamiltonian and thereafter enter the SCF iterations. As a measure of convergence, we adopt the norm of the following matrix in AO basis\n",
"\n",
"$$\n",
"\\mathbf{e} = \n",
"\\mathbf{F D S} -\n",
"\\mathbf{S D F}\n",
"$$\n",
"\n",
"It is convenient to scatter the elements of this matrix into the format of a vector and which we refer to as the *error vector*. During the course of the SCF iterations, we form a sequence of error vectors, $\\mathbf{e}_i$.\n",
"\n",
"The choice of convergence metric may at first appear unintuitive but becomes less so in the MO basis\n",
"\n",
"$$\n",
"\\mathbf{e}^\\mathrm{MO} = \n",
"\\mathbf{C}^\\dagger \\mathbf{e} \\, \\mathbf{C} =\n",
"\\begin{bmatrix}\n",
"0 & -F_\\mathrm{ov} \\\\\n",
"F_\\mathrm{vo} & 0 \\\\\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"where it is seen to correspond to vanishing occupied--virtual blocks in the Fock matrix."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "a6385ce4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter SCF energy Error norm\n",
" 0 -68.84975229 3.09e+00\n",
" 1 -69.95937641 2.88e+00\n",
" 2 -73.34743276 2.83e+00\n",
" 3 -73.46688910 2.23e+00\n",
" 4 -74.74058933 2.18e+00\n",
" 5 -75.55859127 1.41e+00\n",
" 6 -75.86908635 8.26e-01\n",
" 7 -75.97444165 4.82e-01\n",
" 8 -76.00992921 2.74e-01\n",
" 9 -76.02143957 1.57e-01\n",
"10 -76.02519173 8.89e-02\n",
"11 -76.02640379 5.06e-02\n",
"12 -76.02679653 2.88e-02\n",
"13 -76.02692347 1.64e-02\n",
"14 -76.02696455 9.31e-03\n",
"15 -76.02697784 5.30e-03\n",
"16 -76.02698213 3.01e-03\n",
"17 -76.02698352 1.71e-03\n",
"18 -76.02698397 9.74e-04\n",
"19 -76.02698412 5.54e-04\n",
"20 -76.02698416 3.15e-04\n",
"21 -76.02698418 1.79e-04\n",
"22 -76.02698418 1.02e-04\n",
"23 -76.02698419 5.80e-05\n",
"24 -76.02698419 3.30e-05\n",
"25 -76.02698419 1.88e-05\n",
"26 -76.02698419 1.07e-05\n",
"27 -76.02698419 6.07e-06\n",
"28 -76.02698419 3.45e-06\n",
"29 -76.02698419 1.96e-06\n",
"30 -76.02698419 1.12e-06\n",
"31 -76.02698419 6.35e-07\n",
"SCF iterations converged!\n"
]
}
],
"source": [
"max_iter = 50\n",
"conv_thresh = 1e-6\n",
"\n",
"# initial guess from core Hamiltonian\n",
"C = get_MO_coeff(h)\n",
"\n",
"print(\"iter SCF energy Error norm\")\n",
"\n",
"for iter in range(max_iter):\n",
"\n",
" D = np.einsum(\"ik,jk->ij\", C[:, :nocc], C[:, :nocc])\n",
"\n",
" J = np.einsum(\"ijkl,kl->ij\", g, D)\n",
" K = np.einsum(\"ilkj,kl->ij\", g, D)\n",
" F = h + 2 * J - K\n",
"\n",
" E = np.einsum(\"ij,ij->\", h + F, D) + V_nuc\n",
"\n",
" # compute convergence metric\n",
" e_mat = np.linalg.multi_dot([F, D, S]) - np.linalg.multi_dot([S, D, F])\n",
" e_vec = e_mat.reshape(-1)\n",
" error = np.linalg.norm(e_vec)\n",
"\n",
" print(f\"{iter:>2d} {E:16.8f} {error:10.2e}\")\n",
"\n",
" if error < conv_thresh:\n",
" print(\"SCF iterations converged!\")\n",
" break\n",
"\n",
" C = get_MO_coeff(F)"
]
},
{
"cell_type": "markdown",
"id": "cd94e1b9",
"metadata": {},
"source": [
"# Convergence acceleration\n",
"\n",
"## Direct inversion iterative subspace\n",
"\n",
"The Rothaan–Hall scheme suffer from poor numerical convergence and in practice some version of convergence acceleration is adopted. In the method of direct inversion of the iterative subspace (DIIS) {cite}`Pulay1980, Pulay1982, Sellers1993` information is used from not only the present but also previous iterations to form an *averaged* effective one-electron Hamiltonian, or Fock matrix, according to\n",
"\n",
"$$\n",
"\\mathbf{F}_n^\\mathrm{DIIS} = \\sum_{i=1}^n w_i \\mathbf{F}_i\n",
"$$\n",
"\n",
"where $\\mathbf{F}_i$ is the Fock matrix generated in SCF iteration $i$, $n$ is present iteration, and the weights, $w_i$, are to be determined. To guarantee that the one-electron Hamiltonian is preserved in the effective Fock operator, we impose the condition\n",
"\n",
"$$\n",
"\\sum_{i=1}^n w_i = 1\n",
"$$\n",
"\n",
"Under the assumption of a strict linearity between Fock matrices and error vectors, the error vector of the averaged Fock matrix would equal\n",
"\n",
"$$\n",
"\\mathbf{e}_n^\\mathrm{DIIS} = \\sum_{i=1}^n w_i \\mathbf{e}_i\n",
"$$\n",
"\n",
"As the molecular orbitals change from one iteration to the next, this is not strictly the case but it is a good approximation. We can then determine the weights by minimizing the squared norm of this error vector under the imposed constraint. The squared norm becomes equal to\n",
"\n",
"$$\n",
"\\| \\mathbf{e}_n^\\mathrm{DIIS} \\|^2 =\n",
"\\sum_{i,j=1}^n w_i B_{ij} w_j ; \\quad\n",
"B_{ij} = \\langle \\mathbf{e}_i | \\mathbf{e}_j \\rangle\n",
"$$\n",
"\n",
"and the constrained minimization is achieved by introducing a Lagrangian\n",
"\n",
"$$\n",
"L =\n",
"\\| \\mathbf{e}_n^\\mathrm{DIIS} \\|^2 - 2\\lambda\n",
"\\Big(\n",
"\\sum_{i=1}^n w_i - 1\n",
"\\Big)\n",
"$$\n",
"\n",
"where the factor of $-2$ multiplying the Lagrange multiplier $\\lambda$ is a mere convention as to arrive at an explicit matrix equation of the form\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"B_{11} & \\cdots & B_{1n} & -1 \\\\\n",
"\\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
"B_{n1} & \\cdots & B_{nn} & -1 \\\\\n",
"-1 & \\cdots & -1 & 0\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"w_1 \\\\ \\vdots \\\\ w_n \\\\ \\lambda\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"0 \\\\ \\vdots \\\\ 0 \\\\ -1\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"We solve this equation for the weights, $w_i$, and then determine the averaged Fock matrix, $\\mathbf{F}_n^\\mathrm{DIIS}$."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "ad8a7c36",
"metadata": {},
"outputs": [],
"source": [
"def c1diis(F_mats, e_vecs):\n",
"\n",
" n = len(e_vecs)\n",
"\n",
" # build DIIS matrix\n",
" B = -np.ones((n + 1, n + 1))\n",
" B[n, n] = 0\n",
"\n",
" for i in range(n):\n",
" for j in range(n):\n",
" B[i, j] = np.dot(e_vecs[i], e_vecs[j])\n",
"\n",
" b = np.zeros(n + 1)\n",
" b[n] = -1\n",
"\n",
" w = np.matmul(np.linalg.inv(B), b)\n",
"\n",
" F_diis = np.zeros((norb, norb))\n",
" for i in range(n):\n",
" F_diis += w[i] * F_mats[i]\n",
"\n",
" return F_diis"
]
},
{
"cell_type": "markdown",
"id": "12be4eb0",
"metadata": {},
"source": [
"## SCF iterations\n",
"In principle, the only needed modification in the SCF module to implement the DIIS scheme is to replace the original Fock matrix with the weighted averaged counterpart before the determination of the new MO coefficients. But it is also required to save Fock matrices and error vectors from previous SCF iterations. In practice, this extra storage requirement does not severely hamper applications, and in particular so as an optimum stabilitity in the DIIS scheme is experienced with the use of information from a limited number (about 10) of the previous iterations."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "1406fa21",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter SCF energy Error norm\n",
" 0 -68.84975229 3.09e+00\n",
" 1 -69.95937641 2.88e+00\n",
" 2 -75.88501761 8.37e-01\n",
" 3 -75.97016719 5.18e-01\n",
" 4 -76.01483442 2.34e-01\n",
" 5 -76.02676704 3.29e-02\n",
" 6 -76.02698021 3.32e-03\n",
" 7 -76.02698392 8.01e-04\n",
" 8 -76.02698418 2.04e-04\n",
" 9 -76.02698419 6.97e-05\n",
"10 -76.02698419 2.06e-06\n",
"11 -76.02698419 5.43e-07\n",
"SCF iterations converged!\n"
]
}
],
"source": [
"e_vecs = []\n",
"F_mats = []\n",
"\n",
"# initial guess from core Hamiltonian\n",
"C = get_MO_coeff(h)\n",
"\n",
"print(\"iter SCF energy Error norm\")\n",
"\n",
"for iter in range(max_iter):\n",
"\n",
" D = np.einsum(\"ik,jk->ij\", C[:, :nocc], C[:, :nocc])\n",
"\n",
" J = np.einsum(\"ijkl,kl->ij\", g, D)\n",
" K = np.einsum(\"ilkj,kl->ij\", g, D)\n",
" F = h + 2 * J - K\n",
" F_mats.append(F)\n",
"\n",
" E = np.einsum(\"ij,ij->\", h + F, D) + V_nuc\n",
"\n",
" # compute convergence metric\n",
" e_mat = np.linalg.multi_dot([F, D, S]) - np.linalg.multi_dot([S, D, F])\n",
" e_vecs.append(e_mat.reshape(-1))\n",
" error = np.linalg.norm(e_vecs[-1])\n",
"\n",
" print(f\"{iter:>2d} {E:16.8f} {error:10.2e}\")\n",
"\n",
" if error < conv_thresh:\n",
" print(\"SCF iterations converged!\")\n",
" break\n",
"\n",
" F = c1diis(F_mats, e_vecs)\n",
" C = get_MO_coeff(F)"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}