{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "palestinian-medicaid", "metadata": { "tags": [ "remove_cell" ] }, "outputs": [], "source": [ "\"\"\"The Roothaan-Hall self-consistent field procedure\"\"\"\n", "\n", "__author__ = \"Roberto Di Remigio\"\n", "__credit__ = [\"Roberto Di Remigio\", \"Xin Li\"]\n", "\n", "__copyright__ = \"(c) 2021, ENCSS and PDC\"\n", "__license__ = \"MIT\"\n", "__date__ = \"2021-04-13\"" ] }, { "cell_type": "markdown", "id": "natural-flood", "metadata": { "tags": [ "remove_cell" ] }, "source": [ "
\n", " \n", "
" ] }, { "cell_type": "markdown", "id": "solid-thong", "metadata": { "tags": [] }, "source": [ "# The Roothaan-Hall self-consistent field procedure" ] }, { "cell_type": "markdown", "id": "wired-jacksonville", "metadata": { "tags": [ "remove_cell" ] }, "source": [ "
\n", " \n", " \n", "
\n", " To run the selected code cell, hit
Shift + Enter
\n", "
\n", "
" ] }, { "cell_type": "markdown", "id": "robust-selling", "metadata": {}, "source": [ "## Theory refresher\n", "\n", "In Hartree-Fock (HF) theory, we approximate the molecular electronic wavefunction as a single Slater determinant:\n", "\n", "$$\n", "| \\Phi \\rangle = \n", "| \\phi_{1} \\cdots \\phi_{N} \\rangle = \n", " \\frac{1}{\\sqrt{N!}} \\det\n", " \\begin{pmatrix}\n", " \\phi_{1}(x_1) & \\cdots & \\phi_{N}(x_1) \\\\\n", " \\vdots & \\ddots & \\vdots \\\\\n", " \\phi_{1}(x_N) & \\cdots & \\phi_{N}(x_N) \\\\\n", " \\end{pmatrix},\n", "$$\n", "\n", "with $N = N_{\\alpha} + N_{\\beta}$ is the number of electrons and $\\phi_{p}$ the *molecular spin-orbitals* (MOs). From the Slater-Condon rules, the HF energy is:\n", "\n", "$$\n", "E_{\\mathrm{HF}} = \n", "\\sum_{i=1} \\langle \\phi_{i} | -\\frac{1}{2}\\nabla^{2} + V_{\\mathrm{Ne}} | \\phi_{i} \\rangle + \n", "\\frac{1}{2} \\sum_{ij} (\\langle \\phi_{i}\\phi_{j}| \\phi_{i}\\phi_{j} \\rangle - \\langle \\phi_{i}\\phi_{j}| \\phi_{j}\\phi_{i} \\rangle) = \n", "\\sum_{i}h_{ii} + \n", "\\frac{1}{2} \\sum_{ij} \\langle ij | ij \\rangle - \\langle ij | ji \\rangle.\n", "$$" ] }, { "cell_type": "markdown", "id": "stone-disclaimer", "metadata": {}, "source": [ "### Atomic orbital basis expansion\n", "\n", "The MOs will be expanded in a fixed, localized, *nonorthogonal* atomic orbital (AO) basis:\n", "\n", "$$\n", "\\phi^{\\sigma}_{p} = \\sum_{\\mu}^{N_{\\mathrm{AO}}} \\chi_{\\mu}C^{\\sigma}_{\\mu p} \n", "$$\n", "\n", "where $\\sigma$ is the spin of the MO: either $\\alpha$ or $\\beta$. The energy expression becomes:\n", "\n", "$$\n", "E_{\\mathrm{HF}} = \n", "\\sum_{\\mu \\nu} D_{\\mu\\nu} h_{\\mu\\nu} + \n", "\\frac{1}{2}\n", "\\sum_{\\mu\\nu\\kappa\\lambda} \n", "D_{\\mu\\nu}D_{\\kappa\\lambda}(\\mu\\nu|\\kappa\\lambda) - \n", "\\frac{a}{2}\\sum_{\\mu\\nu\\kappa\\lambda} \n", "\\left[ D^{\\alpha}_{\\mu\\kappa}D^{\\alpha}_{\\nu\\lambda} + D^{\\beta}_{\\mu\\kappa}D^{\\beta}_{\\nu\\lambda} \\right] (\\mu\\nu|\\kappa\\lambda)\n", "$$\n", "\n", "with the AO density matrices defined as:\n", "\n", "$$\n", " D^{\\sigma}_{\\mu\\nu} = \\sum_{i=1}^{N_{\\sigma}} C_{\\mu i}^{\\sigma} C_{\\nu i}^{\\sigma}, \\quad D_{\\mu\\nu} = D^{\\alpha}_{\\mu\\nu} + D^{\\beta}_{\\mu\\nu}.\n", "$$\n", "\n", "We will develop our HF code for *closed-shell molecules*, $N_{\\alpha} = N_{\\beta} = \\frac{N}{2}$, and adopt a *spin-restricted formulation*, where $\\alpha$ and $\\beta$ spin-orbitals share the same spatial part:\n", "\n", "$$\n", "D^{\\alpha}_{\\mu\\nu} = D^{\\beta}_{\\mu\\nu} = \\frac{1}{2} D_{\\mu\\nu} = \\sum_{i=1}^{N/2} C_{\\mu i} C_{\\nu i}\n", "$$" ] }, { "cell_type": "markdown", "id": "legitimate-parcel", "metadata": {}, "source": [ "### Variational optimization\n", "\n", "In order to obtain the best possible single-determinant representation in the given AO basis, we optimize the expansion coefficients to obtain a stationary point of the restricted Hatree-Fock energy under the constraint that the MOs remain orthonormal. \n", "We set up a *Lagrangian*:\n", "\n", "$$\n", "L_{\\mathrm{RHF}} = E_{\\mathrm{RHF}} - \\sum_{ij}\\varepsilon_{ij}[\\langle i | j \\rangle - \\delta_{ij} ],\n", "$$\n", "\n", "and seek one of its stationary states:\n", "\n", "$$\n", "\\frac{\\partial L_{\\mathrm{RHF}}}{\\partial C_{\\theta p}} = 0.\n", "$$\n", "\n", "Through some algebra and remembering that the orbitals are *invariant* to unitary transformation *within* the occupied and virtual spaces, this is equivalent to the generalized eigenvalue problem:\n", "\n", "$$\n", "\\mathbf{F}\\mathbf{C} = \\mathbf{S}\\mathbf{C}\\mathbf{\\varepsilon},\n", "$$\n", "\n", "with $\\mathbf{F}$ the *Fock matrix*:\n", "\n", "$$\n", "F_{\\mu\\nu} = \n", "h_{\\mu\\nu} + \n", "\\sum_{\\kappa\\lambda} \\left[ 2(\\mu\\nu|\\kappa\\lambda) - (\\mu\\kappa|\\nu\\lambda) \\right]D_{\\lambda\\kappa} =\n", "h_{\\mu\\nu} + 2J_{\\mu\\nu} - K_{\\mu\\nu}.\n", "$$\n", "\n", "The eigenvalue problem is *nonlinear*: the Fock matrix depends on its eigenvectors through the density matrix." ] }, { "cell_type": "markdown", "id": "contained-vocabulary", "metadata": {}, "source": [ "## The Roothaan-Hall algorithm\n", "\n", "At a glance, the algorithm for solving the **Roothaan-Hall equations** is as follows:\n", "\n", "![Roothaan-Hall self-consistent field iterations](../img/rh-scf.svg)\n", "\n", "We will start our implementation by creating a molecule and assigning it a basis set." ] }, { "cell_type": "code", "execution_count": null, "id": "numerical-glory", "metadata": {}, "outputs": [], "source": [ "import veloxchem as vlx\n", "\n", "h2o_xyz = \"\"\"3\n", "water \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", "mol = vlx.Molecule.from_xyz_string(h2o_xyz)\n", "print(mol.get_string())\n", "\n", "basis = vlx.MolecularBasis.read(mol, \"sto-3g\")\n", "print(basis.get_string(mol))" ] }, { "cell_type": "markdown", "id": "greenhouse-aurora", "metadata": {}, "source": [ "Next, we generate all necessary one-electron integrals: overlap, kinetic, and nuclear potential. Note that we *immediately* convert all the integrals to NumPy arrays using the `to_numpy` function." ] }, { "cell_type": "code", "execution_count": null, "id": "serial-budget", "metadata": {}, "outputs": [], "source": [ "ovldrv = vlx.OverlapIntegralsDriver()\n", "ovl = ovldrv.compute(mol, basis)\n", "S = ovl.to_numpy()" ] }, { "cell_type": "code", "execution_count": null, "id": "handled-tuesday", "metadata": {}, "outputs": [], "source": [ "kindrv = vlx.KineticEnergyIntegralsDriver()\n", "kin = kindrv.compute(mol, basis)\n", "T = kin.to_numpy()" ] }, { "cell_type": "code", "execution_count": null, "id": "effective-flush", "metadata": {}, "outputs": [], "source": [ "naidrv = vlx.NuclearPotentialIntegralsDriver()\n", "nai = naidrv.compute(mol, basis)\n", "V = -nai.to_numpy()" ] }, { "cell_type": "markdown", "id": "thick-cleaning", "metadata": {}, "source": [ "Finally, we generate the ERI 4-index tensor. When using the `compute_in_mem` method, VeloxChem will return a NumPy array directly." ] }, { "cell_type": "code", "execution_count": null, "id": "civilian-brush", "metadata": {}, "outputs": [], "source": [ "eridrv = vlx.ElectronRepulsionIntegralsDriver()\n", "eri = eridrv.compute_in_mem(mol, basis)" ] }, { "cell_type": "markdown", "id": "neutral-notice", "metadata": {}, "source": [ "### Matrix and tensor operations with NumPy\n", "\n", "We will leverage the [NumPy library](https://numpy.org/) for all our linear algebra needs. Before moving on to the implementation of the actual algorithm, we will review the functions that we will use." ] }, { "cell_type": "code", "execution_count": null, "id": "legal-conjunction", "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "id": "choice-stock", "metadata": {}, "source": [ "For matrix multiplications and, in general, tensor contractions, it is very convenient to use [np.einsum](https://numpy.org/devdocs/reference/generated/numpy.einsum.html). This function accepts a string, definining the operation to perform and a variadic pack of operands. The string uses the [Einstein summation convention](https://en.wikipedia.org/wiki/Einstein_notation).\n", "For example, the matrix-matrix multiplication:\n", "\n", "$$\n", "C_{ij} = \\sum_{k}A_{ik}B_{kj}\n", "$$\n", "\n", "is performed with:\n", "\n", "```\n", "C = np.einsum(\"ik,kj->ij\", A, B)\n", "```\n", "\n", "NumPy includes functions for linear algebra, such as decompositions and direct linear solvers. These are in the `linalg` submodule. For example, `np.linalg.eigh` performs the full diagonalization of a Hermitian matrix and return its eigenvalues *and* its eigenvectors:\n", "\n", "```\n", "v, w = np.linalg.eigh(A)\n", "```\n", "\n", "We will also need to take *slices* of matrices and tensors and [advanced indexing](https://numpy.org/doc/stable/reference/arrays.indexing.html) in NumPy will help with this task. For example, to select the columns between indices $i$ and $j$ in a matrix:\n", "\n", "```\n", "slice = U[:, i:j]\n", "```\n", "\n", "the single `:` means to select all rows." ] }, { "cell_type": "markdown", "id": "acquired-subcommittee", "metadata": {}, "source": [ "## Orthogonalization of the AO basis\n", "\n", "The Roothaan-Hall equations are a generalized eigenvalue problem, due to the nonorthogonality of the AO basis. In order to use the `np.linalg.eigh` function from NumPy, we first need to obtain the *orthogonal* AO (OAO) basis, that is, find a transformation matrix $\\mathbf{X}$ such that:\n", "\n", "$$\n", "\\mathbf{X}^{t} \\mathbf{S} \\mathbf{X} = \\mathbf{1}\n", "$$\n", "\n", "In the absence of linear dependencies in the AO basis, the overlap matrix is symmetric (and positive-definite) It is thus diagonalizable in a basis of its eigenvectors:\n", "\n", "$$\n", "\\mathbf{S} = \\mathbf{U}\\mathbf{\\sigma}\\mathbf{U}^{t}\n", "$$\n", "\n", "with $\\mathbf{\\sigma} = \\mathrm{diag}(\\sigma_{0}, \\cdots, \\sigma_{n})$ the diagonal matrix of eigenvalues and $\\mathbf{U}$ its eigenvectors. Any function of the overlap matrix can be defined in terms of its eigenvalues and eigenvectors. This gives us two choices for the orthogonalization:\n", "\n", "- **Symmetric (Löwdin)**, with $\\mathbf{X} = \\mathbf{U}\\mathbf{\\sigma}^{-\\frac{1}{2}}\\mathbf{U}^{t}$.\n", "- **Canonical**, with $\\mathbf{X} = \\mathbf{U}\\mathbf{\\sigma}^{-\\frac{1}{2}}$.\n", "\n", "In both cases, we achieve a factorization of the overlap matrix as $\\mathbf{S} = \\left(\\mathbf{X}^{t}\\mathbf{X}\\right)^{-1}$." ] }, { "cell_type": "code", "execution_count": null, "id": "color-pittsburgh", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "sigma, U = np.linalg.eigh(S)\n", "\n", "# TODO form the X matrix for *symmetric* orthogonalization\n", "X = \n", "\n", "# check that the orthogonalizer is correct\n", "# TODO re-form S from the orthogonalizer matrix\n", "test_S =\n", "np.testing.assert_allclose(test_S, S, atol=1e-10)" ] }, { "cell_type": "markdown", "id": "blind-shirt", "metadata": {}, "source": [ "## Initial guess\n", "\n", "We use the very simple core guess: the initial MO coefficients are the eigenvectors of the core Hamiltonian: $h = T + V$. The algoritm goes as follows:\n", "1. Transform the core Hamiltonian $h$ to OAO basis.\n", "2. Diagonalize it.\n", "3. Backtransform the coefficients to AO basis.\n", "4. Form the density matrix using $D_{\\mu\\nu} = \\sum_{i=1}^{N/2} C_{\\mu i} C_{\\nu i}$." ] }, { "cell_type": "code", "execution_count": null, "id": "willing-messaging", "metadata": {}, "outputs": [], "source": [ "N_O = mol.number_of_electrons() // 2\n", "\n", "# TODO form the core Hamiltonian\n", "h =\n", "# TODO transform to orthogonal AO basis using the orthogonalizer X\n", "h_OAO =\n", "# TODO diagonalize\n", "_, C_OAO =\n", "# TODO backtransform coefficients from OAO basis to AO basis\n", "C = \n", "# TODO form density matrix from occupied slice of coefficients in AO basis\n", "D =" ] }, { "cell_type": "markdown", "id": "allied-assault", "metadata": {}, "source": [ "We can check that the density matrix so produced is idempotent:\n", "\n", "$$\n", "\\mathbf{D}\\mathbf{S}\\mathbf{D} = \\mathbf{D}\n", "$$\n", "\n", "and normalized to the number of electrons:\n", "\n", "$$\n", "N = 2\\,\\mathrm{tr} \\, \\mathbf{D}\\mathbf{S} = 2 \\sum_{\\mu\\nu}D_{\\mu\\nu}S_{\\nu\\mu}\n", "$$\n", "\n", "These checks will ensure that our implementation produces sane density matrices!" ] }, { "cell_type": "code", "execution_count": null, "id": "crucial-colors", "metadata": {}, "outputs": [], "source": [ "# TODO check that the density matrix is normalized to the number of electrons\n", "np.testing.assert_allclose(???, mol.number_of_electrons(), atol=1e-10)\n", "# TODO check that the density matrix is idempotent\n", "np.testing.assert_allclose(???, D, atol=1e-10)" ] }, { "cell_type": "markdown", "id": "behavioral-living", "metadata": {}, "source": [ "The initial energy is:\n", "\n", "$$\n", "E^{[0]} = E_{\\mathrm{NN}} + \\mathrm{tr}\\,\\mathbf{h}\\mathbf{D}^{[0]}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "attended-beatles", "metadata": {}, "outputs": [], "source": [ "E_NN = mol.nuclear_repulsion_energy()\n", "E = np.einsum(\"pq,qp->\", h, D) + E_NN\n", "print(f\"Initial SCF energy: {E:20.12f}\")" ] }, { "cell_type": "markdown", "id": "sophisticated-winter", "metadata": {}, "source": [ "## The iterative procedure\n", "\n", "With a guess density at hand, we can start the Roothaan-Hall self-consistent iterations. We will run until a convergence criterion is met, `conv_thresh`, or until we exceed the maximum number of iterations, `max_iter`.\n", "We choose the *norm of the difference between two successive density matrices*, $\\| \\mathbf{D}^{[n+1]} - \\mathbf{D}^{[n]} \\|$, as convergence criterion." ] }, { "cell_type": "code", "execution_count": null, "id": "english-harvard", "metadata": {}, "outputs": [], "source": [ "max_iter = 50\n", "conv_thresh = 1e-6" ] }, { "cell_type": "markdown", "id": "civilian-spanking", "metadata": {}, "source": [ "At each iteration, we form the *Coulomb*, $J$, and *exchange*, $K$, matrices from a contraction of the density matrix and the pre-computed ERI tensor:\n", "\n", "$$\n", " J^{[n]}_{\\mu\\nu} = \\sum_{\\kappa\\lambda}(\\mu\\nu|\\kappa\\lambda)D_{\\lambda\\kappa}^{[n]}\n", "$$\n", "\n", "and \n", "\n", "$$\n", " K^{[n]}_{\\mu\\nu} = \\sum_{\\kappa\\lambda}(\\mu\\kappa|\\nu\\lambda)D_{\\lambda\\kappa}^{[n]}\n", "$$\n", "\n", "Finally, the energy is:\n", "\n", "$$\n", "E^{[n]} = E_{\\mathrm{NN}} + \n", " \\mathrm{tr}\\,(\\mathbf{h} + \\mathbf{F}^{[n]})\\mathbf{D}^{[n]}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "confused-silence", "metadata": {}, "outputs": [], "source": [ "print(\"# it. SCF energy (E_h) |Delta_E| Delta_D\")\n", "print(f\" 0 {E:20.12f} N/A N/A\")\n", "\n", "for it in range(1, max_iter+1):\n", " # TODO build Fock matrix in AO basis\n", " F = \n", " \n", " # TODO transform Fock matrix to OAO basis using the orthogonalizer X\n", " F_OAO = \n", " \n", " # TODO diagonalize F_OAO\n", " _, C_OAO = \n", " \n", " # TODO backtransform coefficients to AO basis\n", " C =\n", " \n", " # TODO form density matrix from occupied slice of coefficients\n", " D_ =\n", " \n", " # evaluate energy\n", " E_ = np.einsum(\"pq,qp->\", h+F, D_) + E_NN\n", " \n", " # TODO compute convergence metric\n", " Delta_D =\n", " \n", " # print iteration stats\n", " print(f\" {it:>2d} {E_:20.12f} {np.abs(E_ - E):7.5e} {Delta_D:7.5e}\")\n", " \n", " # update current density and current energy\n", " D = D_\n", " E = E_\n", " \n", " # TODO check if convergence threshold was met\n", " \n", " # TODO fail if maximum number of iterations was exceeded without achieving convergence" ] }, { "cell_type": "markdown", "id": "general-truck", "metadata": {}, "source": [ "## Testing our implementation\n", "\n", "Let's compare the result of our own SCF with that obtainable from VeloxChem." ] }, { "cell_type": "code", "execution_count": null, "id": "signal-tribune", "metadata": { "tags": [ "hide_output" ] }, "outputs": [], "source": [ "scfdrv = vlx.ScfRestrictedDriver()\n", "scfdrv.compute(mol, basis)\n", "ref_E = scfdrv.get_scf_energy()" ] }, { "cell_type": "code", "execution_count": null, "id": "administrative-classification", "metadata": {}, "outputs": [], "source": [ "print(f\"Difference with reference result: {E-ref_E:5.8e}\")" ] }, { "cell_type": "markdown", "id": "heated-dynamics", "metadata": {}, "source": [ "## References\n", "\n", "- Lehtola, S.; Blockhuys, F.; Van Alsenoy, C. An Overview of Self-Consistent Field Calculations Within Finite Basis Sets. Molecules 2020, 25 (5). https://doi.org/10.3390/molecules25051218." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.9" } }, "nbformat": 4, "nbformat_minor": 5 }