diff --git a/exercises/ex08_solution/hartree-fock.ipynb b/exercises/ex08_solution/hartree-fock.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..c538bd0183e512b495c89bc4dcb2a7529b0a68fc
--- /dev/null
+++ b/exercises/ex08_solution/hartree-fock.ipynb
@@ -0,0 +1,389 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Hartree-Fock for Electron Ground States of Atoms and Molecules"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import scipy.linalg as la\n",
+ "from numpy import dot, pi, exp, sqrt\n",
+ "from scipy.special import erf\n",
+ "\n",
+ "# interactive plots\n",
+ "#%matplotlib notebook\n",
+ "# nice inline plots\n",
+ "%matplotlib inline\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "plt.rcParams['figure.figsize'] = 16, 9"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## H Atom"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We approximate the single Slater-type orbital with Gaussian functions, using pre-calculated exponents. This is called the STO-nG basis set (with $n=4$ in our case)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "alpha = [13.00773, 1.962079, 0.444529, 0.1219492]\n",
+ "dim = len(alpha)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Set up and solve the generalised eigenvalue problem $\\sum H_{i,j} d_j = \\epsilon \\sum S_{i,j} d_j$ as derived in the script / exercise session."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "def spq(p,q):\n",
+ " return (np.pi/(alpha[p]+alpha[q]))**1.5\n",
+ "\n",
+ "def tpq(p,q):\n",
+ " return 3*alpha[p]*alpha[q]*np.pi**1.5/(alpha[p]+alpha[q])**2.5\n",
+ "\n",
+ "def apq(p,q):\n",
+ " return -2*np.pi/(alpha[p]+alpha[q])\n",
+ "\n",
+ "h = np.array([[tpq(p,q)+apq(p,q) for q in range(dim)] for p in range(dim)]) # Hamiltonian\n",
+ "s = np.array([[spq(p,q) for q in range(dim)] for p in range(dim)]) # overlap matrix\n",
+ "energies = la.eigh(h,s,eigvals_only=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Ground state energy [Hartree]: -0.499278405667\n"
+ ]
+ }
+ ],
+ "source": [
+ "print('Ground state energy [Hartree]:', energies[0])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The ground state for the 1s electron of the H atom is known to be $-0.5 E_h$ ($\\simeq 13.6 eV$), which makes the approximation using four Gaussian-type-orbitals decent."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## He Atom"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Again we refer to a published set of GTO basis function exponential coefficients for the He atom (e.g. from https://bse.pnl.gov)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "alpha = [0.297104, 1.236745, 5.749982, 38.216677]\n",
+ "dim = len(alpha)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Build equation (9) from the exercise sheet, derived as shown in the exercise session."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "def sij(p,q):\n",
+ " \"\"\"overlap matrix elements\"\"\"\n",
+ " return (np.pi/(alpha[p]+alpha[q]))**1.5\n",
+ "\n",
+ "def tij(p,q):\n",
+ " \"\"\"non-interacting matrix elements\"\"\"\n",
+ " return 3*alpha[p]*alpha[q]*np.pi**1.5/(alpha[p]+alpha[q])**2.5 - 4*np.pi/(alpha[p]+alpha[q])\n",
+ "\n",
+ "def vijkl(i,j,k,l):\n",
+ " \"\"\"Hartree matrix elements\"\"\"\n",
+ " return 2*np.pi**2.5/(alpha[i]+alpha[j])/(alpha[k]+alpha[l])/np.sqrt(alpha[i]+alpha[j]+alpha[k]+alpha[l])\n",
+ "\n",
+ "s = np.array([[sij(p,q) for q in range(dim)] for p in range(dim)]) # overlap matrix\n",
+ "t = np.array([[tij(p,q) for q in range(dim)] for p in range(dim)]) # non-interacting matrix\n",
+ "v = np.array([[[[vijkl(i,j,k,l) for i in range(dim)] for j in range(dim)] \n",
+ " for k in range(dim)] for l in range(dim)]) # Hartree matrix\n",
+ "\n",
+ "d = np.ones(dim) # initial coefficient vector\n",
+ "d /= np.sqrt(dot(d,dot(s,d))) # normalize\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note how in the derivation, we have seen terms of $d$ in the Fock matrix $f$.\n",
+ "Since our equation is not a real generalised eigenvalue problem anymore, we need to emulate one through a self-consistency iteration. Convergence depends heavily on the choice of basis, for which STO-4G is known to behave nicely."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "eps = -0.669956328027\n",
+ "eps = -0.669956328027\n",
+ "Ground state energy [Hartree]: -2.07854760879\n"
+ ]
+ }
+ ],
+ "source": [
+ "tol = 1e-10\n",
+ "eps = 0; oldeps = eps+2*tol\n",
+ "while abs(eps-oldeps) > tol:\n",
+ " f = t + dot(dot(v,d),d) # Fock operator matrix\n",
+ " ens,vecs = la.eigh(f,s) # solve GEV problem\n",
+ " oldeps = eps\n",
+ " minidx = np.argmin(ens)\n",
+ " eps = ens[minidx]\n",
+ " d = vecs[:,minidx]\n",
+ " d /= np.sqrt(dot(d,dot(s,d))) # normalize\n",
+ " print('eps =', eps)\n",
+ "\n",
+ "e0 = 2*dot(dot(t,d),d) + dot(dot(dot(dot(v,d),d),d),d)\n",
+ "print('Ground state energy [Hartree]:', e0)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The experimental value is $-2.903$ Hartree. STO-4G already starts to show inadequacies."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## H$_2$ Molecule"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's try to use the exponentials for the GTOs of a single H atom to approximate a H$_2$ molecule."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "alpha = [13.00773, 1.962079, 0.444529, 0.1219492]*2\n",
+ "centr = [0]*4 + [1]*4 # center of basis functions\n",
+ "dim = len(alpha)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Take care to center each set of GTOs around one nucleus at $R_A, R_B$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "alpha = [13.00773, 1.962079, 0.444529, 0.1219492]*2\n",
+ "centr = [0]*4 + [1]*4 # center of basis functions\n",
+ "dim = len(alpha)\n",
+ "\n",
+ "def kexp(i,j):\n",
+ " \"\"\"recurring factor in Gaussian overlap integrals\"\"\"\n",
+ " return exp(-alpha[i]*alpha[j]/(alpha[i]+alpha[j])*(centr[i]-centr[j])**2)\n",
+ "\n",
+ "def rp(i,j):\n",
+ " \"\"\"weighted center position R_P\"\"\"\n",
+ " return (alpha[i]*centr[i] + alpha[j]*centr[j])/(alpha[i]+alpha[j])\n",
+ "\n",
+ "def f0(q):\n",
+ " \"\"\"F_0(q)\"\"\"\n",
+ " if q == 0: return 1\n",
+ " return 0.5*sqrt(pi/q)*erf(sqrt(q))\n",
+ "\n",
+ "def sij(i,j):\n",
+ " \"\"\"overlap matrix elements\"\"\"\n",
+ " return (pi/(alpha[i]+alpha[j]))**1.5 * kexp(i,j)\n",
+ "\n",
+ "def kinij(i,j):\n",
+ " \"\"\"kinetic energy matrix element\"\"\"\n",
+ " a = alpha[i]; b = alpha[j]\n",
+ " return a*b/(a+b) * (3 - 2*a*b/(a+b)*(centr[i]-centr[j])**2) * (pi/(a+b))**1.5 * kexp(i,j)\n",
+ "\n",
+ "def nucij(i,j,rc):\n",
+ " \"\"\"nuclear attraction matrix element for nucleus at position rc\"\"\"\n",
+ " a = alpha[i]; b = alpha[j]\n",
+ " return -2*pi/(a+b) * kexp(i,j) * f0((a+b)*(rp(i,j)-rc)**2)\n",
+ "\n",
+ "def tij(i,j):\n",
+ " \"\"\"non-interacting matrix elements\"\"\"\n",
+ " return kinij(i,j) + nucij(i,j,0) + nucij(i,j,1)\n",
+ "\n",
+ "def vijkl(i,j,k,l):\n",
+ " \"\"\"Hartree matrix elements\"\"\"\n",
+ " aij = alpha[i]+alpha[j]\n",
+ " akl = alpha[k]+alpha[l]\n",
+ " q = aij*akl/(aij+akl) * (rp(i,j) - rp(k,l))**2\n",
+ " return 2*sqrt(aij*akl/pi/(aij+akl))*sij(i,j)*sij(k,l)*f0(q)\n",
+ "\n",
+ "s = np.array([[sij(p,q) for q in range(dim)] for p in range(dim)]) # overlap matrix\n",
+ "t = np.array([[tij(p,q) for q in range(dim)] for p in range(dim)]) # non-interacting matrix\n",
+ "v = np.array([[[[vijkl(i,j,k,l) for i in range(dim)] for j in range(dim)] \n",
+ " for k in range(dim)] for l in range(dim)]) # Hartree matrix\n",
+ "\n",
+ "d = np.ones(dim) # initial coefficient vector\n",
+ "d /= sqrt(dot(d,dot(s,d))) # normalize"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "eps = -0.669956328027\n",
+ "eps = -0.669956328027\n",
+ "Ground state energy [Hartree]: -1.07854760879\n"
+ ]
+ }
+ ],
+ "source": [
+ "tol = 1e-10\n",
+ "eps = 0; oldeps = eps+2*tol\n",
+ "while abs(eps-oldeps) > tol:\n",
+ " f = t + dot(dot(v,d),d) # Fock operator matrix\n",
+ " ens,vecs = la.eigh(f,s) # solve GEV problem\n",
+ " oldeps = eps\n",
+ " minidx = np.argmin(ens)\n",
+ " eps = ens[minidx]\n",
+ " d = vecs[:,minidx]\n",
+ " d /= sqrt(dot(d,dot(s,d))) # normalize\n",
+ " print('eps =', eps)\n",
+ "\n",
+ "e0 = 1 + 2*dot(dot(t,d),d) + dot(dot(dot(dot(v,d),d),d),d)\n",
+ "print('Ground state energy [Hartree]:', e0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Don't forget to add 1 to the resulting energy: during derivation, the repulsion of the H cores (which is normalised to $1$ Hartree) has been left out."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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.4.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/exercises/exercise8.pdf b/exercises/exercise8.pdf
index 822d3cf353deacf12b9f1997c45233cad9825aa9..f1935e489fdcf3d0eab0f575e3f23885a55cb8f5 100644
Binary files a/exercises/exercise8.pdf and b/exercises/exercise8.pdf differ