hartree-fock.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hartree-Fock in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hartee-Fock boils down to a nonlinear generalized eigenvalue equation. You derive it basically by writing down the many-electron Schrodinger equation for a wavefunction that's a single slater determinant formed from molecular orbitals, each of which is expanded in some basis set containing $N$ elements.\n",
    "\n",
    "$$\n",
    "F C = S C e\n",
    "$$\n",
    "\n",
    "$F, S$ and $C$ are $N \\times N$ matrices, where $N$ is the number of basis functions. $e$ is a diagonal matrix containing the orbital energies (eigenvalues).\n",
    "\n",
    "The Fock matrix, $F$ is\n",
    "\n",
    "$$\n",
    "F_{ij} = H_{ij} + \\sum_a \\sum_{kl} C_{ka} C_{al} [2(ij|kl) - (il|kj)] \n",
    "$$\n",
    "\n",
    "This is also written as \n",
    "$$\n",
    "F_{ij} = H_{ij} + 2J - K\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$S$ is the basis function overlap matrix\n",
    "\n",
    "$$\n",
    "S_{ij} = \\int dr \\; \\chi_i(r) \\chi_j(r)\n",
    "$$\n",
    "\n",
    "$H$ is the core hamiltonian, and is built from two terms.\n",
    "\n",
    "$$\n",
    "H_{ij} = T_{ij} + V_{ij}\n",
    "$$\n",
    "\n",
    "where $T$ is the kinetic energy matrix,\n",
    "\n",
    "$$\n",
    "T_{ij} = -\\frac{1}{2} \\int dr\\; \\chi_i(r) \\nabla^2  \\chi_j(r)\n",
    "$$\n",
    "\n",
    "The matrix $V$ comes from the external potential. For the field caused by the nuclei, we have\n",
    "\n",
    "$$\n",
    "V_{ij} =  \\sum_a \\int dr\\; \\chi_i(r) \\frac{Z_a}{|r-r_a|} \\chi_j(r)\n",
    "$$\n",
    "\n",
    "$J$ is the coulomb matrix\n",
    "\n",
    "$$\n",
    "J_{ij} = \\sum_{kl} D_{kl} (ij|kl)\n",
    "$$\n",
    "\n",
    "$K$ is the exchange matrix\n",
    "$$\n",
    "K_{ij} = \\sum_{kl} D_{kl} (il|kj)\n",
    "$$\n",
    "\n",
    "$D$ is the density matrix\n",
    "$$\n",
    "D_{kl} = \\sum_a C_{ka} C_{al}\n",
    "$$\n",
    "\n",
    "The two-electron integrals are a 4-index tensor\n",
    "\n",
    "$$\n",
    "(ij|kl) = \\int dr_1 dr_2 \\; \\chi_i(r_1) \\chi_j(r_1) \\frac{1}{|r_1-r_2|} \\chi_k(r_2) \\chi_l(r_2)\n",
    "$$\n",
    "\n",
    "For quantum chemistry, the most common types of basis functions are Gaussian of the form\n",
    "$$\n",
    "\\chi(\\mathbf{r}) = (r_x - c_x)^{a_x} (r_y - c_y)^{a_y} (r_z - c_z)^{a_z} \\exp(-\\alpha \\cdot |\\mathbf{r} - \\mathbf{c}|^2)\n",
    "$$\n",
    "\n",
    "The primary reason that these basis functions are widely used is that you can compute the two-electron integrals very efficiently. The equations are quite messy though."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rmcgibbo/miniconda/envs/2.7/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    }
   ],
   "source": [
    "# pip install git+https://github.com/rpmuller/pyquante2.git\n",
    "from __future__ import print_function, division, absolute_import\n",
    "from pyquante2.basis.basisset import basisset\n",
    "from pyquante2.ints.integrals import onee_integrals, twoe_integrals\n",
    "\n",
    "import numpy as np\n",
    "from scipy.linalg import eigh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def rhf(molecule, basis, n_iters=50):\n",
    "    energies = []\n",
    "    bfs = basisset(molecule, basis)\n",
    "    onee = onee_integrals(bfs, molecule)\n",
    "    twoe_tensor = twoe_integrals(bfs)._2e_ints\n",
    "\n",
    "    H = onee.T + onee.V\n",
    "    S = onee.S\n",
    "    e, C = eigh(H, S)\n",
    "\n",
    "    D = np.dot(C[:, :molecule.nclosed()], C[:, :molecule.nclosed()].T)\n",
    "\n",
    "    for i in range(n_iters):\n",
    "        J = np.einsum('kl,ijkl->ij', D, twoe_tensor)\n",
    "        K = np.einsum('kl,ilkj->ij', D, twoe_tensor)\n",
    "        F = H + 2*J - K\n",
    "\n",
    "        e, C = eigh(F, S)\n",
    "        D = np.dot(C[:, :molecule.nclosed()], C[:, :molecule.nclosed()].T)\n",
    "\n",
    "        energies.append(molecule.nuclear_repulsion() + np.sum((2*H + 2*J-K) * D))\n",
    "\n",
    "    return energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stoichiometry = H2O, Charge = 0, Multiplicity = 1\n",
      "8 O     0.000000     0.000000     0.091686\n",
      "1 H     1.422968     0.000000    -0.981210\n",
      "1 H    -1.422968     0.000000    -0.981210\n",
      "Energy -75.9843038988\n"
     ]
    }
   ],
   "source": [
    "from pyquante2.geo.samples import h2o\n",
    "print(h2o)\n",
    "energies = rhf(h2o, '6-31g')\n",
    "print('Energy', energies[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# http://cccbdb.nist.gov/energy3x.asp?method=1&basis=9&charge=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}