added to ex08 solution: drawbacks

9bfab327 · Donjan Rodic · 3a110410 · 9bfab327
Commit 9bfab327 authored 8 years ago by Donjan Rodic
--- a/exercises/ex08_solution/hartree-fock.ipynb
+++ b/exercises/ex08_solution/hartree-fock.ipynb
@@ -352,7 +352,29 @@
   "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."
+    "Don't forget to add 1 to the resulting energy: during derivation, the repulsion of the H cores ( $\\frac{1}{|R_A-R_B|}=1$ Hartree in the Hamiltonian) has been dropped."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "## Drawbacks"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "* The Born-Oppenheimer approximation neglects nuclei, although this doesn't factor too strongly into the examples we have used.\n",
+    "\n",
+    "* The choice of a too small basis and/or one consisting of too simple functions (here GTOs) can be insufficient to describe the system accurately. Bigger sets (STO-6G) and/or more elaborate functions get computationally expensive quickly.\n",
+    "\n",
+    "* The wave function is approximated by a single Slater determinant. Each electron sees the average density of all other electrons (compare: mean field theory), which doesn't factor in electron correlation.\n",
+    "\n",
+    "The last point is addressed extensively in \"post-Hartree-Fock\" methods (e.g. using a sum of \"excited\" Slater determinants)."
   ]
  },
  {

 %% Cell type:markdown id: tags:
 # Hartree-Fock for Electron Ground States of Atoms and Molecules
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 import scipy.linalg as la
 from numpy import dot, pi, exp, sqrt
 from scipy.special import erf
 # interactive plots
 #%matplotlib notebook
 # nice inline plots
 %matplotlib inline
 import matplotlib.pyplot as plt
 plt.rcParams['figure.figsize'] = 16, 9
 ```
 %% Cell type:markdown id: tags:
 ## H Atom
 %% Cell type:markdown id: tags:
 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 id: tags:
 ``` python
 alpha = [13.00773, 1.962079, 0.444529, 0.1219492]
 dim = len(alpha)
 ```
 %% Cell type:markdown id: tags:
 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 id: tags:
 ``` python
 def spq(p,q):
    return (np.pi/(alpha[p]+alpha[q]))**1.5
 def tpq(p,q):
    return 3*alpha[p]*alpha[q]*np.pi**1.5/(alpha[p]+alpha[q])**2.5
 def apq(p,q):
    return -2*np.pi/(alpha[p]+alpha[q])
 h = np.array([[tpq(p,q)+apq(p,q) for q in range(dim)] for p in range(dim)]) # Hamiltonian
 s = np.array([[spq(p,q)          for q in range(dim)] for p in range(dim)]) # overlap matrix
 energies = la.eigh(h,s,eigvals_only=True)
 ```
 %% Cell type:code id: tags:
 ``` python
 print('Ground state energy [Hartree]:', energies[0])
 ```
 %% Output
    Ground state energy [Hartree]: -0.499278405667
 %% Cell type:markdown id: tags:
 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 id: tags:
 ## He Atom
 %% Cell type:markdown id: tags:
 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 id: tags:
 ``` python
 alpha = [0.297104, 1.236745, 5.749982, 38.216677]
 dim = len(alpha)
 ```
 %% Cell type:markdown id: tags:
 Build equation (9) from the exercise sheet, derived as shown in the exercise session.
 %% Cell type:code id: tags:
 ``` python
 def sij(p,q):
    """overlap matrix elements"""
    return (np.pi/(alpha[p]+alpha[q]))**1.5
 def tij(p,q):
    """non-interacting matrix elements"""
    return 3*alpha[p]*alpha[q]*np.pi**1.5/(alpha[p]+alpha[q])**2.5 - 4*np.pi/(alpha[p]+alpha[q])
 def vijkl(i,j,k,l):
    """Hartree matrix elements"""
    return 2*np.pi**2.5/(alpha[i]+alpha[j])/(alpha[k]+alpha[l])/np.sqrt(alpha[i]+alpha[j]+alpha[k]+alpha[l])
 s = np.array([[sij(p,q) for q in range(dim)] for p in range(dim)]) # overlap matrix
 t = np.array([[tij(p,q) for q in range(dim)] for p in range(dim)]) # non-interacting matrix
 v = np.array([[[[vijkl(i,j,k,l) for i in range(dim)] for j in range(dim)]
                                for k in range(dim)] for l in range(dim)]) # Hartree matrix
 d = np.ones(dim) # initial coefficient vector
 d /= np.sqrt(dot(d,dot(s,d))) # normalize
 ```
 %% Cell type:markdown id: tags:
 Note how in the derivation, we have seen terms of $d$ in the Fock matrix $f$.
 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 id: tags:
 ``` python
 tol = 1e-10
 eps = 0; oldeps = eps+2*tol
 while abs(eps-oldeps) > tol:
    f = t + dot(dot(v,d),d) # Fock operator matrix
    ens,vecs = la.eigh(f,s) # solve GEV problem
    oldeps = eps
    minidx = np.argmin(ens)
    eps = ens[minidx]
    d  = vecs[:,minidx]
    d /= np.sqrt(dot(d,dot(s,d))) # normalize
    print('eps =', eps)
 e0 = 2*dot(dot(t,d),d) + dot(dot(dot(dot(v,d),d),d),d)
 print('Ground state energy [Hartree]:', e0)
 ```
 %% Output
    eps = -0.669956328027
    eps = -0.669956328027
    Ground state energy [Hartree]: -2.07854760879
 %% Cell type:markdown id: tags:
 The experimental value is $-2.903$ Hartree. STO-4G already starts to show inadequacies.
 %% Cell type:markdown id: tags:
 ## H$_2$ Molecule
 %% Cell type:markdown id: tags:
 Let's try to use the exponentials for the GTOs of a single H atom to approximate a H$_2$ molecule.
 %% Cell type:code id: tags:
 ``` python
 alpha = [13.00773, 1.962079, 0.444529, 0.1219492]*2
 centr = [0]*4 + [1]*4   # center of basis functions
 dim = len(alpha)
 ```
 %% Cell type:markdown id: tags:
 Take care to center each set of GTOs around one nucleus at $R_A, R_B$.
 %% Cell type:code id: tags:
 ``` python
 alpha = [13.00773, 1.962079, 0.444529, 0.1219492]*2
 centr = [0]*4 + [1]*4   # center of basis functions
 dim = len(alpha)
 def kexp(i,j):
    """recurring factor in Gaussian overlap integrals"""
    return exp(-alpha[i]*alpha[j]/(alpha[i]+alpha[j])*(centr[i]-centr[j])**2)
 def rp(i,j):
    """weighted center position R_P"""
    return (alpha[i]*centr[i] + alpha[j]*centr[j])/(alpha[i]+alpha[j])
 def f0(q):
    """F_0(q)"""
    if q == 0: return 1
    return 0.5*sqrt(pi/q)*erf(sqrt(q))
 def sij(i,j):
    """overlap matrix elements"""
    return (pi/(alpha[i]+alpha[j]))**1.5 * kexp(i,j)
 def kinij(i,j):
    """kinetic energy matrix element"""
    a = alpha[i]; b = alpha[j]
    return a*b/(a+b) * (3 - 2*a*b/(a+b)*(centr[i]-centr[j])**2) * (pi/(a+b))**1.5 * kexp(i,j)
 def nucij(i,j,rc):
    """nuclear attraction matrix element for nucleus at position rc"""
    a = alpha[i]; b = alpha[j]
    return -2*pi/(a+b) * kexp(i,j) * f0((a+b)*(rp(i,j)-rc)**2)
 def tij(i,j):
    """non-interacting matrix elements"""
    return kinij(i,j) + nucij(i,j,0) + nucij(i,j,1)
 def vijkl(i,j,k,l):
    """Hartree matrix elements"""
    aij = alpha[i]+alpha[j]
    akl = alpha[k]+alpha[l]
    q = aij*akl/(aij+akl) * (rp(i,j) - rp(k,l))**2
    return 2*sqrt(aij*akl/pi/(aij+akl))*sij(i,j)*sij(k,l)*f0(q)
 s = np.array([[sij(p,q) for q in range(dim)] for p in range(dim)]) # overlap matrix
 t = np.array([[tij(p,q) for q in range(dim)] for p in range(dim)]) # non-interacting matrix
 v = np.array([[[[vijkl(i,j,k,l) for i in range(dim)] for j in range(dim)]
                                for k in range(dim)] for l in range(dim)]) # Hartree matrix
 d = np.ones(dim) # initial coefficient vector
 d /= sqrt(dot(d,dot(s,d))) # normalize
 ```
 %% Cell type:code id: tags:
 ``` python
 tol = 1e-10
 eps = 0; oldeps = eps+2*tol
 while abs(eps-oldeps) > tol:
    f = t + dot(dot(v,d),d) # Fock operator matrix
    ens,vecs = la.eigh(f,s) # solve GEV problem
    oldeps = eps
    minidx = np.argmin(ens)
    eps = ens[minidx]
    d  = vecs[:,minidx]
    d /= sqrt(dot(d,dot(s,d))) # normalize
    print('eps =', eps)
 e0 = 1 + 2*dot(dot(t,d),d) + dot(dot(dot(dot(v,d),d),d),d)
 print('Ground state energy [Hartree]:', e0)
 ```
 %% Output
    eps = -0.669956328027
    eps = -0.669956328027
    Ground state energy [Hartree]: -1.07854760879
 %% Cell type:markdown id: tags:
-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.
+Don't forget to add 1 to the resulting energy: during derivation, the repulsion of the H cores ( $\frac{1}{|R_A-R_B|}=1$ Hartree in the Hamiltonian) has been dropped.
+%% Cell type:markdown id: tags:
+## Drawbacks
+%% Cell type:markdown id: tags:
+* The Born-Oppenheimer approximation neglects nuclei, although this doesn't factor too strongly into the examples we have used.
+* The choice of a too small basis and/or one consisting of too simple functions (here GTOs) can be insufficient to describe the system accurately. Bigger sets (STO-6G) and/or more elaborate functions get computationally expensive quickly.
+* The wave function is approximated by a single Slater determinant. Each electron sees the average density of all other electrons (compare: mean field theory), which doesn't factor in electron correlation.
+The last point is addressed extensively in "post-Hartree-Fock" methods (e.g. using a sum of "excited" Slater determinants).
 %% Cell type:code id: tags:
 ``` python
 ```