Tutorials/05_Moller-Plesset/5a_conventional-mp2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Tutorial: Second-Order Moller--Plesset Perturbation Theory (MP2)\"\"\"\n",
    "\n",
    "__author__    = \"Dominic A. Sirianni\"\n",
    "__credit__    = [\"Dominic A. Sirianni\", \"Daniel G. A. Smith\"]\n",
    "\n",
    "__copyright__ = \"(c) 2014-2018, The Psi4NumPy Developers\"\n",
    "__license__   = \"BSD-3-Clause\"\n",
    "__date__      = \"2017-05-23\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Second-Order Moller-Plesset Perturbation Theory (MP2)\n",
    "\n",
    "Moller-Plesset perturbation theory [also referred to as many-body perturbation theory (MBPT)] is an adaptation of the more general Rayleigh-Schrodinger perturbation theory (RSPT), applied to problems in molecular electronic structure theory.  This tutorial will provide a brief overview of both RSPT and MBPT, before walking through an implementation of second-order Moller-Plesset perturbation theory (specifically referred to as MP2) which uses conventional, 4-index ERIs.  \n",
    "\n",
    "### I. Overview of Rayleigh-Schrodinger Perturbation Theory\n",
    "Given the Hamiltonian operator $\\hat{H}$ for a system, perturbation theory solves the Schrodinger equation for that system by rewriting $\\hat{H}$ as\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{H} = \\hat{H}{}^{(0)} + \\lambda\\hat{V},\n",
    "\\tag{[Szabo:1996], pp. 322, Eqn. 6.3}\n",
    "\\end{equation}\n",
    "\n",
    "were $\\hat{H}{}^{(0)}$ is the Hamiltonian operator corresponding to a solved problem which resembles $\\hat{H}$, and $\\hat{V}$ is the *perturbation* operator, defined as $\\hat{V} = \\hat{H} - \\hat{H}{}^{(0)}$.  Then the Schrodinger equation for the system becomes \n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{H}\\mid\\Psi_n\\rangle = (\\hat{H}{}^{(0)} + \\lambda\\hat{V})\\mid\\Psi_n\\rangle = E_n\\mid\\Psi_n\\rangle.\n",
    "\\tag{[Szabo:1996], pp. 322, Eqn. 6.2}\n",
    "\\end{equation}\n",
    "\n",
    "The energies $E_n$ and wavefunctions $\\mid\\Psi_n\\rangle$ will both be functions of $\\lambda$; they can therefore be written as a Taylor series expansion about $\\lambda = 0$ ([Szabo:1996], pp. 322, Eqns. 6.4a & 6.4b):\n",
    "\n",
    "\\begin{align}\n",
    "E_n &= E_n^{(0)} + \\lambda E_n^{(1)} + \\lambda^2E_n^{(2)} + \\ldots;\\tag{[Szabo:1996], pp. 322, Eqn. 6.4a}\\\\\n",
    "\\mid\\Psi_n\\rangle &= \\mid\\Psi_n^{(0)}\\rangle + \\lambda\\mid\\Psi_n^{(1)}\\rangle + \\lambda^2\\mid\\Psi_n^{(2)}\\rangle + \\ldots,\\tag{[Szabo:1996], pp. 322, Eqn. 6.4b}\\\\\n",
    "\\end{align}\n",
    "\n",
    "in practice, these perturbation expansions may be truncated to a given power of $\\lambda$.  Substituting the perturbation series above back into the Schrodinger equation yields\n",
    "\n",
    "\\begin{equation*}\n",
    "(\\hat{H}{}^{(0)} + \\lambda\\hat{V})(\\mid\\Psi_n^{(0)}\\rangle + \\lambda\\mid\\Psi_n^{(1)}\\rangle + \\lambda^2\\mid\\Psi_n^{(2)}\\rangle + \\ldots) = (E_n^{(0)} + \\lambda E_n^{(1)} + \\lambda^2E_n^{(2)} + \\ldots)(\\mid\\Psi_n^{(0)}\\rangle + \\lambda\\mid\\Psi_n^{(1)}\\rangle + \\lambda^2\\mid\\Psi_n^{(2)}\\rangle + \\ldots),\n",
    "\\end{equation*}\n",
    "\n",
    "which by equating powers of $\\lambda$ ([Szabo:1996], pp. 323, Eqns. 6.7a-6.7d) gives expressions for the $E_n^{(i)}$and $\\mid\\Psi_n^{(i)}\\rangle$. Note that for $\\lambda^0$, $E_n^{(0)}$ and $\\mid\\Psi_n^{(0)}\\rangle$ are known, as they are the solution to the zeroth-order problem $\\hat{H}{}^{(0)}$.  For $\\lambda^1$ and $\\lambda^2$, the expressions for $E_n^{(1)}$ and $E_n^{(2)}$ are given by\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda^1:\\;\\;\\;\\;E_n^{(1)} &= \\langle\\Psi_n^{(0)}\\mid\\hat{V}\\mid\\Psi_n^{(0)}\\rangle,\\tag{[Szabo:1996], pp. 323, Eqn. 6.8b}\\\\\n",
    "\\lambda^2:\\;\\;\\;\\;E_n^{(2)} &= \\sum_{\\mu\\neq n}\\frac{\\mid\\langle\\Psi_{n}^{(0)}\\mid\\hat{V}\\mid\\Psi_{\\mu}^{(0)}\\rangle\\mid^2}{E_n^{(0)} - E_{\\mu}^{(0)}}\\tag{[Szabo:1996], pp. 324, Eqn. 6.12}\\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II. Overview of Moller-Plesset Perturbation Theory\n",
    "\n",
    "The exact electronic Hamiltonian for an N-electron molecule with a given atomic configuration is (in atomic units):\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{H}_{elec} = \\sum_i\\hat{h}(i) + \\sum_{i<j}\\frac{1}{r_{ij}}.\\tag{[Szabo:1996], pp. 43, Eqn. 2.10}\n",
    "\\end{equation}\n",
    "\n",
    "Moller-Plesset perturbation theory seeks to solve the molecular electronic Schrodinger equation with the above Hamiltonian using the techniques of Rayleigh-Schroding perturbation theory, by selecting the zeroth-order Hamiltonian and wavefunctions to be those from Hartree-Fock theory:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{H}{}^{(0)}  &= \\sum_i \\hat{f}(i) = \\sum_i \\hat{h}(i) + \\upsilon^{HF}(i)\\tag{[Szabo:1996], pp. 350, Eqn. 6.59}\\\\\n",
    "\\hat{V} &= \\hat{H} - \\hat{H}{}^{(0)} = \\sum_{i<j} \\frac{1}{r_{ij}} - \\upsilon^{HF}(i).\\tag{[Szabo:1996], pp. 350, Eqn. 6.60}\\\\\n",
    "\\end{align}\n",
    "\n",
    "With these choices of $\\hat{H}{}^{(0)}$, $\\hat{V}$, and $\\mid\\Psi_n^{(0)}\\rangle$, the ground-state zeroth-order energy is given by $E_0^{(0)} = \\sum_i\\epsilon_i$ ([Szabo:1996], pp. 351, Eqn. 6.67).  Then, the first-order ground state energy may be computed by $E_0^{(1)} = \\langle\\Psi_0^{HF}\\mid\\hat{V}\\mid\\Psi_0^{HF}\\rangle$ to find that the total ground-state energy through first order, $E_0 = E_0^{(0)} + E_0^{(1)}$, is exactly the Hartree-Fock energy. ([Szabo:1996], pp. 351, Eqn. 6.69)  Therefore, the first correction to the Hartree-Fock energy will occur at second-order in the perturbation series, i.e., with $E_0^{(2)}$; truncating the perturbation series at second order is commonly referred to as MP2.  The second order correction to the ground state energy will be given by\n",
    "\n",
    "\\begin{equation}\n",
    "E_0^{(2)} = \\sum_{\\mu\\neq 0}\\frac{\\left|\\langle\\Psi_{0}^{HF}\\mid\\hat{V}\\mid\\Psi_{\\mu}^{HF}\\rangle\\right|^2}{E_0^{(0)} - E_{\\mu}^{(0)}}\\tag{[Szabo:1996], pp. 351, Eqn. 6.70}\n",
    "\\end{equation}\n",
    "\n",
    "For brevity, we will now drop the \"HF\" from all zeroth-order wavefunctions.  This summation is over the eigenstate index $\\mu$, each associated with a different eigenstate of the zeroth-order Hamiltonian.  For a single-determinant wavefunction constructed from spin orbitals, the summation over the eigenstate index $\\mu\\neq 0$ therefore refers to determinants which are constructed from *different* spin orbitals than the ground state determinant.  To distinguish such determinants, we will denote MOs occupied in the ground state with indices $i,\\,j,\\,\\ldots$, and MOs which are unoccupied in the ground state with indices $a,\\,b,\\,\\ldots\\,$.  Then a determinant where orbital $a$ is substituted for orbital $i$ is denoted $\\mid\\Psi_i^a\\rangle$, and so on.  Before substituting this new notation into the above energy expression, however, we may immediately recognize that many terms $\\langle\\Psi_{0}\\mid\\hat{V}\\mid\\Psi_{\\mu}\\rangle$ will not contribute to the second order energy:\n",
    "\n",
    "| Term         | Determinant                               |  Contribution to $E_0^{(2)}$                            |\n",
    "|--------------|-------------------------------------------|---------------------------------------------------------|\n",
    "| Singles      | $\\mid\\Psi_i^a\\rangle$                     | 0; Brillouin's Theorem                                  |\n",
    "| Doubles      | $\\mid\\Psi_{ij}^{ab}\\rangle$               | Survive                                                 |\n",
    "| Higher-order | $\\mid\\Psi_{ijk\\ldots}^{abc\\ldots}\\rangle$ | 0; $\\hat{V}$ is a two-particle operator            |\n",
    "\n",
    "Hence we see that only doubly-substituted determinants $\\mid\\Psi_{ij}^{ab}\\rangle$ will contribute to $E_0^{(2)}$.  From Hartree-Fock theory, we know that \n",
    "\n",
    "\\begin{equation}\n",
    "\\langle\\Psi_0\\mid r_{ij}^{-1}\\mid\\Psi_{ij}^{ab}\\rangle = [ia\\| jb],\\tag{[Szabo:1996], pp. 72, Tbl. 2.6}\n",
    "\\end{equation}\n",
    "\n",
    "where $[ia\\| jb] = [ia\\mid jb] - [ib\\mid ja]$ is an antisymmetrized two-electron integral, and the square brackets \"$[\\ldots ]$\" indicate that we are employing chemists' notation.  What about the energies of these doubly substituted determinants?  Recognizing that the difference between the energies of the newly- and formerly-occupied orbitals in each substitution must modify the total energy of the ground state determinant, \n",
    "\n",
    "\\begin{equation}\n",
    "E_{ij}^{ab} = E_0 - (\\epsilon_i - \\epsilon_a + \\epsilon_j - \\epsilon_b).\n",
    "\\end{equation}\n",
    "\n",
    "Substituting these expressions into the one for the second-order energy, we have that\n",
    "\n",
    "\\begin{equation}\n",
    "E_0^{(2)} = \\sum_{i<j}\\sum_{a<b} \\frac{\\left|\\,[ia\\|jb]\\,\\right|^2}{\\epsilon_i - \\epsilon_a + \\epsilon_j - \\epsilon_b}.\\tag{[Szabo:1996], pp. 352, Eqn. 6.71}\n",
    "\\end{equation}\n",
    "\n",
    "So far, our discussion has used spin-orbitals instead of the more familiar spatial orbitals.  Indeed, significant speedups are achieved when using spatial orbitals.  Integrating out the spin variable $\\omega$ from $E_0^{(2)}$ yields two expressions; one each for the interaction of particles with the same spin (SS) and opposite spin (OS):\n",
    "\n",
    "\\begin{align}\n",
    "E_{\\rm 0,\\,SS}^{(2)} = \\sum_{ij}\\sum_{ab}\\frac{(ia\\mid jb)[(ia\\mid jb) - (ib\\mid ja)]}{\\epsilon_i - \\epsilon_a + \\epsilon_j - \\epsilon_b},\n",
    "E_{\\rm 0,\\,OS}^{(2)} = \\sum_{ij}\\sum_{ab}\\frac{(ia\\mid jb)(ia\\mid jb)}{\\epsilon_i - \\epsilon_a + \\epsilon_j - \\epsilon_b},\n",
    "\\end{align}\n",
    "\n",
    "where these spin-free expressions make use of integrals in chemists' notation over spatial orbitals. (Rearranged from [Szabo:1996], pp. 352, Eqn. 6.74)  Note that an exchange integral arises between particles of the same spin; this is because the motions of particles with identical spins are correlated due to the requirement that $\\left|\\Psi\\right|^2$ remain invariant to the exchange of the spatial and spin coordinates of any pair of electrons.  Finally, the total MP2 correction energy $E_0^{(2)} = E_{\\rm 0,\\,SS}^{(2)} + E_{\\rm 0,\\,OS}^{(2)}$.\n",
    "\n",
    "### Implementation of Conventional MP2\n",
    "\n",
    "Let's begin by importing Psi4 and NumPy, and setting memory and output file options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==> Import statements & Global Options <==\n",
    "import psi4\n",
    "import numpy as np\n",
    "\n",
    "psi4.set_memory(int(2e9))\n",
    "numpy_memory = 2\n",
    "psi4.core.set_output_file('output.dat', False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we can define our molecule and Psi4 options.  Notice that we are using `scf_type pk` to indicate that we wish to use conventional, full 4-index ERIs, and that we have specified `mp2_type conv` so that the MP2 algorithm we check against also uses the conventional ERIs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==> Molecule & Psi4 Options Definitions <==\n",
    "mol = psi4.geometry(\"\"\"\n",
    "O\n",
    "H 1 1.1\n",
    "H 1 1.1 2 104\n",
    "symmetry c1\n",
    "\"\"\")\n",
    "\n",
    "\n",
    "psi4.set_options({'basis':        '6-31g',\n",
    "                  'scf_type':     'pk',\n",
    "                  'mp2_type':     'conv',\n",
    "                  'e_convergence': 1e-8,\n",
    "                  'd_convergence': 1e-8})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since MP2 is a perturbation on the zeroth-order Hartree-Fock description of a molecular system, all of the relevant information (Fock matrix, orbitals, orbital energies) about the system can be computed using any Hartree-Fock program.  We could use the RHF program that we wrote in tutorial 3a, but we could just as easily use Psi4 to do our dirty work.  In the cell below, use Psi4 to compute the RHF energy and wavefunction, and store them using the `return_wfn=True` keyword argument to `psi4.energy()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the SCF wavefunction & energies\n",
    "scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the expression for $E_0^{(2)}$, the two summations are over occupied and virtual indices, respectively.  Therefore, we'll need to get the number of occupied orbitals and the total number of orbitals.  Additionally, we must obtain the MO energy eigenvalues; again since the sums are over occupied and virtual orbitals, it is good to separate the occupied orbital energies from the virtual orbital energies.  From the SCF wavefunction you generated above, get the number of doubly occupied orbitals, number of molecular orbitals, and MO energies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==> Get orbital information & energy eigenvalues <==\n",
    "# Number of Occupied orbitals & MOs\n",
    "ndocc = scf_wfn.nalpha()\n",
    "nmo = scf_wfn.nmo()\n",
    "\n",
    "# Get orbital energies, cast into NumPy array, and separate occupied & virtual\n",
    "eps = np.asarray(scf_wfn.epsilon_a())\n",
    "e_ij = eps[:ndocc]\n",
    "e_ab = eps[ndocc:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the orbital information, Psi4 does not return the ERIs when it does a computation.  Fortunately, however, we can just build them again using the `psi4.core.MintsHelper()` class.  Recall that these integrals will be generated in the AO basis; before using them in the $E_0^{(2)}$ expression, we must transform them into the MO basis.  To do this, we first need to obtain the orbital coefficient matrix, **C**.  In the cell below, generate the ERIs for our molecule, get **C** from the SCF wavefunction, and obtain occupied- and virtual-orbital slices of **C** for future use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Size of the ERI tensor will be 0.00 GB.\n"
     ]
    }
   ],
   "source": [
    "# ==> ERIs <==\n",
    "# Create instance of MintsHelper class\n",
    "mints = psi4.core.MintsHelper(scf_wfn.basisset())\n",
    "\n",
    "# Memory check for ERI tensor\n",
    "I_size = (nmo**4) * 8.e-9\n",
    "print('\\nSize of the ERI tensor will be %4.2f GB.' % I_size)\n",
    "memory_footprint = I_size * 1.5\n",
    "if I_size > numpy_memory:\n",
    "    psi4.core.clean()\n",
    "    raise Exception(\"Estimated memory utilization (%4.2f GB) exceeds allotted memory \\\n",
    "                     limit of %4.2f GB.\" % (memory_footprint, numpy_memory))\n",
    "\n",
    "# Build ERI Tensor\n",
    "I = np.asarray(mints.ao_eri())\n",
    "\n",
    "# Get MO coefficients from SCF wavefunction\n",
    "C = np.asarray(scf_wfn.Ca())\n",
    "Cocc = C[:, :ndocc]\n",
    "Cvirt = C[:, ndocc:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to transform the four-index integrals from the AO to the MO basis, we must perform the following contraction:\n",
    "\n",
    "$$(i\\,a\\mid j\\,b) = C_{\\mu i}C_{\\nu a}(\\mu\\,\\nu\\mid|\\,\\lambda\\,\\sigma)C_{\\lambda j}C_{\\sigma b}$$\n",
    "\n",
    "Again, here we are using $i,\\,j$ as occupied orbital indices and $a,\\, b$ as virtual orbital indices.  We could carry out the above contraction all in one step using either `np.einsum()` or explicit loops:\n",
    "\n",
    "~~~python\n",
    "# Naive Algorithm for ERI Transformation\n",
    "Imo = np.einsum('pi,qa,pqrs,rj,sb->iajb', Cocc, Cvirt, I, Cocc, Cvirt, optimize=True)\n",
    "~~~\n",
    "\n",
    "Notice that the transformation from AO index to occupied (virtual) MO index requires only the occupied (virtual) block of the **C** matrix; this allows for computational savings in large basis sets, where the virtual space can be very large.  This algorithm, while simple, has the distinct disadvantage of horrendous scaling.  Examining the contraction more closely, we see that there are 8 unique indices, and thus the step above scales as ${\\cal O}(N^8)$.  With this algorithm, a twofold increase of the number of MO's would result in $2^8 = 256\\times$ expense to perform.  We can, however, refactor the above contraction such that\n",
    "\n",
    "$$(i\\,a\\mid j\\,b) = \\left[C_{\\mu i}\\left[C_{\\nu a}\\left[C_{\\lambda j}\\left[C_{\\sigma b}(\\mu\\,\\nu\\mid|\\,\\lambda\\,\\sigma)\\right]\\right]\\right]\\right],$$\n",
    "\n",
    "where we have now written the transfomation as four ${\\cal O}(N^5)$ steps instead of one ${\\cal O}(N^8)$ step. This is a savings of $\\frac{4}{n^3}$, and is responsible for the feasibility of the MP2 method for application to any but very small systems and/or basis sets.  We may carry out the above ${\\cal O}(N^5)$ algorithm by carrying out one index transformation at a time, and storing the result in a temporary array.  In the cell below, transform the ERIs from the AO to MO basis, using our smarter algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==> Transform I -> I_mo @ O(N^5) <==\n",
    "tmp = np.einsum('pi,pqrs->iqrs', Cocc, I, optimize=True)\n",
    "tmp = np.einsum('qa,iqrs->iars', Cvirt, tmp, optimize=True)\n",
    "tmp = np.einsum('iars,rj->iajs', tmp, Cocc, optimize=True)\n",
    "I_mo = np.einsum('iajs,sb->iajb', tmp, Cvirt, optimize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note here that we can use infrastructure in Psi4 to carry out the above integral transformation; this entails obtaining the occupied and virtual blocks of **C** Psi4-side, and then using the built-in `MintsHelper` function `MintsHelper.mo_eri()` to transform the integrals.  Just to check your work above, execute the next cell to see this tech in action:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Do our transformed ERIs match Psi4's? True\n"
     ]
    }
   ],
   "source": [
    "# ==> Compare our Imo to MintsHelper <==\n",
    "Co = scf_wfn.Ca_subset('AO','OCC')\n",
    "Cv = scf_wfn.Ca_subset('AO','VIR')\n",
    "MO = np.asarray(mints.mo_eri(Co, Cv, Co, Cv))\n",
    "print(\"Do our transformed ERIs match Psi4's? %s\" % np.allclose(I_mo, np.asarray(MO)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have all the pieces needed to compute $E_0^{(2)}$.  This could be done by writing explicit loops over occupied and virtual indices Python side, e.g.,\n",
    "~~~python\n",
    "# Compute SS MP2 Correlation\n",
    "mp2_ss_corr = 0.0\n",
    "for i in xrange(ndocc):\n",
    "    for a in xrange(nmo - ndocc):\n",
    "        for j in xrange(ndocc):\n",
    "            for b in xrange(nmo - ndocc):\n",
    "                numerator = I_mo[i,a,j,b] * (I_mo[i, a, j, b] - I_mo[i, b, j, a])\n",
    "                mp2_ss_corr += numerator / (e_ij[i] + e_ij[j] - e_ab[a] - e_ab[b])\n",
    "~~~\n",
    "This method, while very clear what is going on and easy to program, has the distinct disadvantage that Python loops are much slower than the same block written in a compiled language like C, C++, or Fortran.  For this reason, it is better to use infrastructure available in NumPy like `np.einsum()`, `np.dot`, etc. to explicitly compute the above quantity C-side.  It should be clear how to contract the four-index integrals $(i\\,a\\mid j\\,b)$ and $(i\\,a\\mid j\\,b)$ with one another, but what about the energy eigenvalues $\\epsilon$?  We can use a NumPy trick called *broadcasting* to construct a four-index array of all possible energy denominators, which can then be contracted with the full I_mo arrays.  To do this, we'll use the function `np.reshape()`:\n",
    "~~~python\n",
    "# Prepare 4d energy denominator array\n",
    "e_denom = e_ij.reshape(-1, 1, 1, 1) # Diagonal of 4d array are occupied orbital energies\n",
    "e_denom -= e_ab.reshape(-1, 1, 1)   # all combinations of (e_ij - e_ab)\n",
    "e_denom += e_ij.reshape(-1, 1)      # all combinations of [(e_ij - e_ab) + e_ij]\n",
    "e_denom -= e_ab                     # All combinations of full denominator\n",
    "e_denom = 1 / e_denom               # Take reciprocal for contracting with numerator\n",
    "~~~\n",
    "In the cell below, compute the energy denominator using `np.reshape()` and contract this array with the four-index ERIs to compute the same-spin and opposite-spin MP2 correction using `np.einsum()`.  Then, add these quantities to the SCF energy computed above to obtain the total MP2 energy.\n",
    "\n",
    "Hint: For the opposite-spin correlation, use `np.swapaxes()` to obtain the correct ordering of the indices in the exchange integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==> Compute MP2 Correlation & MP2 Energy <==\n",
    "# Compute energy denominator array\n",
    "e_denom = 1 / (e_ij.reshape(-1, 1, 1, 1) - e_ab.reshape(-1, 1, 1) + e_ij.reshape(-1, 1) - e_ab)\n",
    "\n",
    "# Compute SS & OS MP2 Correlation with Einsum\n",
    "mp2_os_corr = np.einsum('iajb,iajb,iajb->', I_mo, I_mo, e_denom, optimize=True)\n",
    "mp2_ss_corr = np.einsum('iajb,iajb,iajb->', I_mo, I_mo - I_mo.swapaxes(1,3), e_denom, optimize=True)\n",
    "\n",
    "# Total MP2 Energy\n",
    "MP2_E = scf_e + mp2_os_corr + mp2_ss_corr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    MP2 Energy........................................................PASSED\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# ==> Compare to Psi4 <==\n",
    "psi4.compare_values(psi4.energy('mp2'), MP2_E, 6, 'MP2 Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## References\n",
    "\n",
    "1. Original paper: \"Note on an Approximation Treatment for Many-Electron Systems\"\n",
    "\t> [[Moller:1934:618](https://journals.aps.org/pr/abstract/10.1103/PhysRev.46.618)] C. Møller and M. S. Plesset, *Phys. Rev.* **46**, 618 (1934)\n",
    "2. The Laplace-transformation in MP theory: \"Minimax approximation for the decomposition of energy denominators in Laplace-transformed Møller–Plesset perturbation theories\"\n",
    "    > [[Takasuka:2008:044112](http://aip.scitation.org/doi/10.1063/1.2958921)] A. Takatsuka, T. Siichiro, and W. Hackbusch, *J. Phys. Chem.*, **129**, 044112 (2008)\n",
    "3. Equations taken from:\n",
    "\t> [[Szabo:1996](https://books.google.com/books?id=KQ3DAgAAQBAJ&printsec=frontcover&dq=szabo+%26+ostlund&hl=en&sa=X&ved=0ahUKEwiYhv6A8YjUAhXLSCYKHdH5AJ4Q6AEIJjAA#v=onepage&q=szabo%20%26%20ostlund&f=false)] A. Szabo and N. S. Ostlund, *Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory*. Courier Corporation, 1996.\n",
    "4. Algorithms taken from:\n",
    "\t> [Crawford:prog] T. D. Crawford, \"The Second-Order Møller–Plesset Perturbation Theory (MP2) Energy.\"  Accessed via the web at http://github.com/CrawfordGroup/ProgrammingProjects."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}