Coupled Cluster with Singles and Doubles (CCSD) in Python

Coupled cluster methods are among the most accurate electronic structure methods available today. For example, with a good choice for a basis, the CCSD(T) equations will give you the correct energy of a molecular system to chemical accuracy (~0.1 kcal/mol). They are also quite formidable to derive, and (in my opinion) to program. Which is why I only coded up CCSD!

Thankfully there are some wonderful resources available to understand the coupled cluster methods. I highly recommend Crawford and Schaefer’s “ An Introduction to Coupled Cluster Theory for Computational Chemists” (2000). I have yet to find a clearer and more complete explanation of coupled cluster theory. The derivations inside it are nasty, but once you get a grasp of diagrammatic techniques, it isn’t so bad :)

In order to understand coupled cluster a bit better, I recently programmed the CCSD energy and amplitude equations in Python. It is for a HeH+ molecule with a bond length of 0.9295 Angstrom, using an STO-3G basis — same system I’ve used before on this blog. The results match what Gaussian09 calculates as well, so I was pretty happy to see it work. As always, I’ve hard-coded the two-electron integrals and other SCF results into the program, so you can just focus on what CCSD does. The functions will look esoteric, and unless you’ve worked with coupled-cluster before, the program should NOT look intuitive or easy to understand — point is, don’t panic. But I’ve provided a reference Stanton (1991) that contains all the equations used. Between Stanton and Crawford, you can understand what is going on here. Read through the comments to get a better idea: the main idea is to take the results of an SCF calculation and apply them to a similarity transformation of the Hamiltonian. The transformed Hamiltonian now contains ‘excited’ determinants, which is the requirement for electron correlation — in other words you get a multi-reference quality calculation from a single reference (a Hartree-Fock calculation).

Anyway, here you have it:

(Special thanks to Junhao Li for pointing out bugs in the code! Everything should be correct now.)

#!/usr/bin/python

####################################
#
#  CCSD ENERGY CALCULATION ON 
#    HeH+ / STO-3G / R = 0.9295 Ang
#
#  Reference for Equations:
#    Stanton, et al
#      J. Chem. Phys 94 (6),
#      15 March 1991
#
####################################

from __future__ import division
import math
import numpy as np

####################################
#
#   FUNCTIONS
#
####################################

# Return compound index given four indices
def eint(a,b,c,d):
  if a > b: ab = a*(a+1)/2 + b
  else: ab = b*(b+1)/2 + a
  if c > d: cd = c*(c+1)/2 + d
  else: cd = d*(d+1)/2 + c
  if ab > cd: abcd = ab*(ab+1)/2 + cd
  else: abcd = cd*(cd+1)/2 + ab
  return abcd

# Return Value of spatial MO two electron integral
# Example: (12\vert 34) = tei(1,2,3,4)
def teimo(a,b,c,d):
  return ttmo.get(eint(a,b,c,d),0.0e0)

####################################
#
#  INITIALIZE ORBITAL ENERGIES 
#  AND TRANSFORMED TWO ELECTRON
#  INTEGRALS  
#
####################################

Nelec = 2 # we have 2 electrons in HeH+
dim = 2 # we have two spatial basis functions in STO-3G
E = [-1.52378656, -0.26763148] # molecular orbital energies
# python dictionary containing two-electron repulsion integrals
ttmo = {5.0: 0.94542695583037617, 12.0: 0.17535895381500544, 14.0: 0.12682234020148653, 17.0: 0.59855327701641903, 19.0: -0.056821143621433257, 20.0: 0.74715464784363106}
ENUC = 1.1386276671 # nuclear repulsion energy for HeH+ -- constant
EN   = -3.99300007772 # SCF energy
####################################################
#
#  CONVERT SPATIAL TO SPIN ORBITAL MO
#
####################################################

# This makes the spin basis double bar integral (physicists' notation)

spinints=np.zeros((dim*2,dim*2,dim*2,dim*2))
for p in range(1,dim*2+1):
  for q in range(1,dim*2+1):
    for r in range(1,dim*2+1):
      for s in range(1,dim*2+1):
        value1 = teimo((p+1)//2,(r+1)//2,(q+1)//2,(s+1)//2) * (p%2 == r%2) * (q%2 == s%2)
        value2 = teimo((p+1)//2,(s+1)//2,(q+1)//2,(r+1)//2) * (p%2 == s%2) * (q%2 == r%2)
        spinints[p-1,q-1,r-1,s-1] = value1 - value2

#####################################################
#
#  Spin basis fock matrix eigenvalues 
#
#####################################################

fs = np.zeros((dim*2))
for i in range(0,dim*2):
    fs[i] = E[i//2]
fs = np.diag(fs) # put MO energies in diagonal array

#######################################################
#
#   CCSD CALCULATION
#
#######################################################
dim = dim*2 # twice the dimension of spatial orbital

# Init empty T1 (ts) and T2 (td) arrays

ts = np.zeros((dim,dim))
td = np.zeros((dim,dim,dim,dim))

# Initial guess T2 --- from MP2 calculation!

for a in range(Nelec,dim):
  for b in range(Nelec,dim):
    for i in range(0,Nelec):
      for j in range(0,Nelec):
        td[a,b,i,j] += spinints[i,j,a,b]/(fs[i,i] + fs[j,j] - fs[a,a] - fs[b,b])

# Make denominator arrays Dai, Dabij
# Equation (12) of Stanton
Dai = np.zeros((dim,dim))
for a in range(Nelec,dim):
  for i in range(0,Nelec):
    Dai[a,i] = fs[i,i] - fs[a,a]

# Stanton eq (13)
Dabij = np.zeros((dim,dim,dim,dim))
for a in range(Nelec,dim):
  for b in range(Nelec,dim):
    for i in range(0,Nelec):
      for j in range(0,Nelec):
        Dabij[a,b,i,j] = fs[i,i] + fs[j,j] - fs[a,a] - fs[b,b]

# Stanton eq (9)
def taus(a,b,i,j):
  taus = td[a,b,i,j] + 0.5*(ts[a,i]*ts[b,j] - ts[b,i]*ts[a,j])
  return taus

# Stanton eq (10)
def tau(a,b,i,j):
  tau = td[a,b,i,j] + ts[a,i]*ts[b,j] - ts[b,i]*ts[a,j]
  return tau

# We need to update our intermediates at the beginning, and 
# at the end of each iteration. Each iteration provides a new
# guess at the amplitudes T1 (ts) and T2 (td), that *hopefully*
# converges to a stable, ground-state, solution.

def updateintermediates(x):
  if x == True:
    # Stanton eq (3)
    Fae = np.zeros((dim,dim))
    for a in range(Nelec,dim):
      for e in range(Nelec,dim):
        Fae[a,e] = (1 - (a == e))*fs[a,e]
        for m in range(0,Nelec):
          Fae[a,e] += -0.5*fs[m,e]*ts[a,m]
          for f in range(Nelec,dim):
            Fae[a,e] += ts[f,m]*spinints[m,a,f,e] 
            for n in range(0,Nelec):
              Fae[a,e] += -0.5*taus(a,f,m,n)*spinints[m,n,e,f]

    # Stanton eq (4)
    Fmi = np.zeros((dim,dim))
    for m in range(0,Nelec):
      for i in range(0,Nelec):
        Fmi[m,i] = (1 - (m == i))*fs[m,i]
        for e in range(Nelec,dim):
          Fmi[m,i] += 0.5*ts[e,i]*fs[m,e]
          for n in range(0,Nelec):
            Fmi[m,i] += ts[e,n]*spinints[m,n,i,e] 
            for f in range(Nelec,dim):
              Fmi[m,i] += 0.5*taus(e,f,i,n)*spinints[m,n,e,f]

    # Stanton eq (5)
    Fme = np.zeros((dim,dim))
    for m in range(0,Nelec):
      for e in range(Nelec,dim):
        Fme[m,e] = fs[m,e]
        for n in range(0,Nelec):
          for f in range(Nelec,dim):
            Fme[m,e] += ts[f,n]*spinints[m,n,e,f]

    # Stanton eq (6)
    Wmnij = np.zeros((dim,dim,dim,dim))
    for m in range(0,Nelec):
      for n in range(0,Nelec):
        for i in range(0,Nelec):
          for j in range(0,Nelec):
            Wmnij[m,n,i,j] = spinints[m,n,i,j]
            for e in range(Nelec,dim):
              Wmnij[m,n,i,j] += ts[e,j]*spinints[m,n,i,e] - ts[e,i]*spinints[m,n,j,e]
              for f in range(Nelec,dim):
                Wmnij[m,n,i,j] += 0.25*tau(e,f,i,j)*spinints[m,n,e,f]

    # Stanton eq (7)
    Wabef = np.zeros((dim,dim,dim,dim))
    for a in range(Nelec,dim):
      for b in range(Nelec,dim):
        for e in range(Nelec,dim):
          for f in range(Nelec,dim):
            Wabef[a,b,e,f] = spinints[a,b,e,f]
            for m in range(0,Nelec):
              Wabef[a,b,e,f] += -ts[b,m]*spinints[a,m,e,f] + ts[a,m]*spinints[b,m,e,f]
              for n in range(0,Nelec):
                Wabef[a,b,e,f] += 0.25*tau(a,b,m,n)*spinints[m,n,e,f]

    # Stanton eq (8)
    Wmbej = np.zeros((dim,dim,dim,dim))
    for m in range(0,Nelec):
      for b in range(Nelec,dim):
        for e in range(Nelec,dim):
          for j in range(0,Nelec):
            Wmbej[m,b,e,j] = spinints[m,b,e,j]
            for f in range(Nelec,dim):
              Wmbej[m,b,e,j] += ts[f,j]*spinints[m,b,e,f]
            for n in range(0,Nelec):
              Wmbej[m,b,e,j] += -ts[b,n]*spinints[m,n,e,j]
              for f in range(Nelec,dim):
                Wmbej[m,b,e,j] += -(0.5*td[f,b,j,n] + ts[f,j]*ts[b,n])*spinints[m,n,e,f]

    return Fae, Fmi, Fme, Wmnij, Wabef, Wmbej

# makeT1 and makeT2, as they imply, construct the actual amplitudes necessary for computing
# the CCSD energy (or computing an EOM-CCSD Hamiltonian, etc)

# Stanton eq (1)
def makeT1(x,ts,td):
  if x == True:
    tsnew = np.zeros((dim,dim))
    for a in range(Nelec,dim):
      for i in range(0,Nelec):
        tsnew[a,i] = fs[i,a]
        for e in range(Nelec,dim):
          tsnew[a,i] += ts[e,i]*Fae[a,e]
        for m in range(0,Nelec):
          tsnew[a,i] += -ts[a,m]*Fmi[m,i]
          for e in range(Nelec,dim):
            tsnew[a,i] += td[a,e,i,m]*Fme[m,e]
            for f in range(Nelec,dim):
              tsnew[a,i] += -0.5*td[e,f,i,m]*spinints[m,a,e,f]
            for n in range(0,Nelec):
              tsnew[a,i] += -0.5*td[a,e,m,n]*spinints[n,m,e,i]
        for n in range(0,Nelec):
          for f in range(Nelec,dim): 
            tsnew[a,i] += -ts[f,n]*spinints[n,a,i,f]
        tsnew[a,i] = tsnew[a,i]/Dai[a,i]
  return tsnew

# Stanton eq (2)
def makeT2(x,ts,td):
  if x == True:
    tdnew = np.zeros((dim,dim,dim,dim))
    for a in range(Nelec,dim):
      for b in range(Nelec,dim):
        for i in range(0,Nelec):
          for j in range(0,Nelec):
            tdnew[a,b,i,j] += spinints[i,j,a,b]
            for e in range(Nelec,dim):
              tdnew[a,b,i,j] += td[a,e,i,j]*Fae[b,e] - td[b,e,i,j]*Fae[a,e]
              for m in range(0,Nelec):
                tdnew[a,b,i,j] += -0.5*td[a,e,i,j]*ts[b,m]*Fme[m,e] + 0.5*td[b,e,i,j]*ts[a,m]*Fme[m,e]
                continue
            for m in range(0,Nelec):
              tdnew[a,b,i,j] += -td[a,b,i,m]*Fmi[m,j] + td[a,b,j,m]*Fmi[m,i]
              for e in range(Nelec,dim):
                tdnew[a,b,i,j] += -0.5*td[a,b,i,m]*ts[e,j]*Fme[m,e] + 0.5*td[a,b,j,m]*ts[e,i]*Fme[m,e]
                continue
            for e in range(Nelec,dim):
              tdnew[a,b,i,j] += ts[e,i]*spinints[a,b,e,j] - ts[e,j]*spinints[a,b,e,i]
              for f in range(Nelec,dim):
                tdnew[a,b,i,j] += 0.5*tau(e,f,i,j)*Wabef[a,b,e,f]
                continue
            for m in range(0,Nelec):
              tdnew[a,b,i,j] += -ts[a,m]*spinints[m,b,i,j] + ts[b,m]*spinints[m,a,i,j]  
              for e in range(Nelec,dim):
                tdnew[a,b,i,j] +=  td[a,e,i,m]*Wmbej[m,b,e,j] - ts[e,i]*ts[a,m]*spinints[m,b,e,j]
                tdnew[a,b,i,j] += -td[a,e,j,m]*Wmbej[m,b,e,i] + ts[e,j]*ts[a,m]*spinints[m,b,e,i]
                tdnew[a,b,i,j] += -td[b,e,i,m]*Wmbej[m,a,e,j] + ts[e,i]*ts[b,m]*spinints[m,a,e,j]
                tdnew[a,b,i,j] +=  td[b,e,j,m]*Wmbej[m,a,e,i] - ts[e,j]*ts[b,m]*spinints[m,a,e,i]
                continue
              for n in range(0,Nelec):
                tdnew[a,b,i,j] += 0.5*tau(a,b,m,n)*Wmnij[m,n,i,j]
                continue
            tdnew[a,b,i,j] = tdnew[a,b,i,j]/Dabij[a,b,i,j] 
    return tdnew

# Expression from Crawford, Schaefer (2000) 
# DOI: 10.1002/9780470125915.ch2
# Equation (134) and (173)
# computes CCSD energy given T1 and T2
def ccsdenergy():
  ECCSD = 0.0
  for i in range(0,Nelec):
    for a in range(Nelec,dim):
      ECCSD += fs[i,a]*ts[a,i]
      for j in range(0,Nelec):
        for b in range(Nelec,dim):
          ECCSD += 0.25*spinints[i,j,a,b]*td[a,b,i,j] + 0.5*spinints[i,j,a,b]*(ts[a,i])*(ts[b,j]) 
  return ECCSD

#================
# MAIN LOOP
# CCSD iteration
#================
ECCSD = 0
DECC = 1.0
while DECC > 0.000000001: # arbitrary convergence criteria
  OLDCC = ECCSD
  Fae,Fmi,Fme,Wmnij,Wabef,Wmbej = updateintermediates(True)
  tsnew = makeT1(True,ts,td)
  tdnew = makeT2(True,ts,td)
  ts = tsnew
  td = tdnew
  ECCSD = ccsdenergy()
  DECC = abs(ECCSD - OLDCC)

print "E(corr,CCSD) = ", ECCSD
print "E(CCSD) = ", ECCSD + ENUC + EN

Boom: -2.8626 Hartrees :)

Equation of Motion (EOM) derivation of RPA

The Equation of Motion derivations of excited state and response properties are elegant, and (in my opinion) very direct. Here I’ll give an example of how it can be used to derive TDHF/RPA. We’ve already done it two other ways, and the agreement between them is important if we wish to extend these methods further.

Given an exact ground state, we can say that

Define operator \({Q_n^{\dagger}}\) and \({Q_n}\):
\(\vert n \rangle = Q_n^{\dagger} \vert 0 \rangle, \qquad Q_n\vert 0 \rangle = 0 \implies Q_n^{\dagger} = \vert n \rangle \langle 0 \vert \ \ \ \ \\)

These operators generate excited states from the ground state (not excited determinants, as in the case of post-HF correlation methods). So it is clear that, when acting on an exact ground state:

Multiply on left by arbitrary state of form \({\langle 0 \vert \delta Q}\), giving
\(\langle 0 \vert [\delta Q, [H, Q_n^{\dagger}]]\vert 0 \rangle = \hbar\omega_{0n} \langle 0 \vert [\delta Q, Q_n^{\dagger}] \vert 0 \rangle \ \ \ \ \\)

Where we have made use of the fact that \({\langle 0 \vert Q_n^{\dagger} = \langle 0 \vert HQ_n^{\dagger} = 0}\). Note that is we express \({Q}\) by particle-hole operators \({a_p^{\dagger}a_q}\), \({a_p^{\dagger}a_q^{\dagger}a_r a_s}\) with coefficients \({C_{pq}}\) and \({C_{pqrs}}\), then \({\delta Q}\) is given by \({\frac{\partial Q}{\partial C} \delta C}\) for arbitrary variations \({\delta C}\). These are in principle exact, since \({\delta Q \vert 0 \rangle}\) exhausts the whole Hilbert space, such that the above equation corresponds to the full Schrödinger equation. Tamm-Dancoff (or Configuration Interaction Singles) can be obtained by approximating \({\vert 0 \rangle \rightarrow \vert \hbox{HF} \rangle}\) and the operator \({Q_n^{\dagger} = \sum\limits_{ia} C_{ia} a_a^{\dagger}a_i}\), restricting ourselves to 1p-1h excitations. Thus \({\delta Q \vert 0 \rangle = \sum_{ia} a_a^{\dagger}a_i \vert \hbox{HF} \rangle \delta C_{ai}}\), (\({\delta C_{ai}}\) cancels), and
\(\sum\limits_{bj} \langle \hbox{HF} \vert [a_i^{\dagger}a_a, [H, a_b^{\dagger}a_j]]\vert \hbox{HF} \rangle C_{jb} = \hbar\omega \langle \hbox{HF} \vert [a_i^{\dagger}a_a, a_a^{\dagger}a_i] \vert \hbox{HF} \rangle C_{ia} \ \ \ \ \\)

These are the CIS equations. Put another way:
\(\sum\limits_{bj} \left\{(\epsilon_a - \epsilon_i)\delta_{ab}\delta_{ij} + \langle aj \vert \vert ib \rangle \right\} C_{bj} = E^{CIS}C_{ai} \ \ \ \ \\)

Similarly, for RPA/TDHF, if we consider a ground state containing 2p-2h correlations, we can not only create a p-h pair, but also destroy one. Thus (choosing the minus sign for convenience):
\(Q_n^{\dagger} = \sum\limits_{ia} X_{ia} a_a^{\dagger}a_i - \sum\limits_{ia} Y_{ia} a_i^{\dagger}a_a, \qquad \hbox{and} Q_n \vert RPA \rangle = 0 \ \ \ \ \\)

So instead of the basis of only single excitations, and therefore one matrix \({C_{ia}}\), we work in a basis of single excitations and single de-excitations, and have two matrices \({X_{ia}}\) and \({Y_{ia}}\). We also have two kinds of variations \({\delta Q \vert 0 \rangle}\), namely \({a_a^{\dagger}a_i \vert 0 \rangle}\) and \({a_i^{\dagger}a_a \vert 0 \rangle}\). This gives us two sets of equations:
\(\langle \hbox{RPA} \vert [a_i^{\dagger}a_a, [H,Q_n^{\dagger}]]\vert \hbox{RPA} \rangle = \hbar \omega \langle \hbox{RPA} \vert [a_i^{\dagger}a_a, Q_n^{\dagger}] \vert \hbox{RPA} \rangle \notag \ \ \ \ \\)

These contain only expectation values of our four Fermion operators, which cannot be calculated since we still do not know \({\vert \hbox{RPA} \rangle}\). Thus we assume \({\vert \hbox{RPA} \rangle \rightarrow \vert \hbox{HF} \rangle}\). This gives
\(\langle \hbox{RPA} \vert [a_i^{\dagger}a_a, a_b^{\dagger}a_j] \vert \hbox{RPA} \rangle = \langle \hbox{HF} \vert [a_i^{\dagger}a_a, a_b^{\dagger}a_j] \vert \hbox{HF} \rangle = \delta_{ij}\delta_{ab} \ \ \ \ \\)

The probability of finding states \({a_a^{\dagger}a_i \vert 0 \rangle}\) and \({a_i^{\dagger}a_a \vert 0 \rangle}\) in excited state \({\vert n\rangle}\), that is, the p-h and h-p matrix elements of transition density matrix \({\rho^{(1)}}\) are:
\(\rho^{(1)}_{ai} = \langle 0 \vert a_i^{\dagger}a_a \vert n \rangle \simeq \langle \hbox{HF} \vert [a_i^{\dagger}a_a, Q_n^{\dagger}] \vert \hbox{HF} \rangle = X_{ia} \ \ \ \ \\)

Thus altogether
\(\begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \\ \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \\ \end{bmatrix} = \hbar\omega \begin{bmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{-1} \\ \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \\ \end{bmatrix} \ \ \ \ \\)

with
\(A_{ia,jb} = \langle \hbox{HF} \vert [a_i^{\dagger}a_a, [H,a_b^{\dagger}a_j]]\vert \hbox{HF} \rangle \ \ \ \ \\)

and
\(B_{ia,jb} = - \langle \hbox{HF} \vert [a_i^{\dagger}a_a, [H,a_j^{\dagger}a_b]]\vert \hbox{HF} \rangle \ \ \ \ \\)

which are the TDHF/RPA equations.

Polarization Propagator Derivation of TDHF/RPA

Following immediately where we left off last time, when \({A}\) and \({B}\) are one-electron number conserving, e.g. creation/annihilation operators are in pairs, we have

which allow us to define the general polarization propagator. If the elements did not conserve the number of electrons, we could use the same formalism to get the electron propagator, which can describe the attachment/detachment of electrons from a system, thus describing processes such as ionization. In general, it is sufficient to consider the typical elements \({A \rightarrow a_r^{\dagger}a_s}\) and \({B \rightarrow a_t^{\dagger}a_u}\). If spin-orbitals of reference are \({\psi_i, \psi_j}\),…, then a particle-hole basis includes only excitation and de-excitation operators \({a_a^{\dagger}a_i}\) and \({a_i^{\dagger}a_a}\), respectively. Thus the basis and its dual:
\(\mathbf{n} = (\mathbf{e} \quad \mathbf{d}), \qquad \mathbf{n}^{\dagger} = (\mathbf{e} \quad \mathbf{d})^{\dagger} \ \ \ \ \\)
\(e_{ia} = a_a^{\dagger}a_i = E_{ia}, \qquad d_{ai} = a_i^{\dagger}a_a = E_{ia}^{\dagger} \ \ \ \ \ (29)\)

This gives, for our resolvent,
\(R(\omega) = \left[\mathbf{n}^{\dagger}(\hbar\omega\hat{1} - \hat{H})\mathbf{n}\right]^{-1} = \begin{bmatrix} \mathbf{e}^{\dagger}(\hbar\omega\hat{1} - \hat{H})\mathbf{e} & \mathbf{e}^{\dagger}(\hbar\omega\hat{1} - \hat{H})\mathbf{d} \\ \mathbf{d}^{\dagger}(\hbar\omega\hat{1} - \hat{H})\mathbf{e} & \mathbf{d}^{\dagger}(\hbar\omega\hat{1} - \hat{H})\mathbf{d} \\ \end{bmatrix}^{-1} \ \ \ \ \ (30)\)

As an example, we give the elements of the first block:
\(E^{\dagger}_{ia}(\hbar\omega\hat{1} - \hat{H})E_{jb} = \hbar\omega E^{\dagger}_{ia} E_{jb} - E^{\dagger}_{ia} \left[H,E_{jb}\right] \ \ \ \ \ (31)\)

with scalar products
\(E^{\dagger}_{ia} E_{jb} = \langle \psi_0 \vert [E^{\dagger}_{ia},E_{jb}]\vert \psi_0 \rangle, \qquad E^{\dagger}_{ia} \left[H,E_{jb}\right] = \langle \psi_0 \vert [E^{\dagger}_{ia},[H,E_{jb}]]\vert \psi_0 \rangle \ \ \ \ \ (32)\)

Evaluating these products gives
\(E^{\dagger}_{ia} E_{jb} = \langle \psi_0 \vert [E^{\dagger}_{ia},E_{jb}]\vert \psi_0 \rangle = \langle \psi_0 \vert E^{\dagger}_{ia}E_{jb}\vert \psi_0 \rangle - \langle \psi_0 \vert E_{jb}E^{\dagger}_{ia}\vert \psi_0 \rangle = \langle \psi_0 \vert E^{\dagger}_{ia}E_{jb}\vert \psi_0 \rangle = \delta_{ia,jb} \ \ \ \ \ (33)\)

and

\(E^{\dagger}_{ia} \left[H,E_{jb}\right] = \langle \psi_0 \vert [E^{\dagger}_{ia},[H,E_{jb}]]\vert \psi_0 \rangle \ \ \ \ \\)
\(E^{\dagger}_{ia} \left[H,E_{jb}\right] = \langle \psi_0 \vert E^{\dagger}_{ia} H E_{jb}\vert \psi_0 \rangle - \langle \psi_0 \vert E^{\dagger}_{ia}E_{jb} H\vert \psi_0 \rangle - \langle \psi_0 \vert H E_{jb} E^{\dagger}_{ia}\vert \psi_0 \rangle + \langle \psi_0 \vert E_{jb} H E^{\dagger}_{ia} \vert \psi_0 \rangle \ \ \ \ \\)
\(E^{\dagger}_{ia} \left[H,E_{jb}\right] = \langle \psi^a_i \vert H \vert \psi_j^b \rangle - \delta_{ia}\delta_{jb}\langle \psi_0 \vert H \vert \psi_0 \rangle \ \ \ \ \\)
\(E^{\dagger}_{ia} \left[H,E_{jb}\right] = E_0\delta_{ij}\delta_{ab} + F_{ab}\delta_{ij} - F_{ij}\delta_{ab} + \langle aj \vert \vert ib \rangle - \delta_{ia}\delta{jb}E_0 \ \ \ \ \\)
\(E^{\dagger}_{ia} \left[H,E_{jb}\right] = (\epsilon_a - \epsilon_i)\delta_{ia}\delta_{ab} + \langle aj \vert \vert ib \rangle = \mathbf{A} \ \ \ \ \\)
Which is precisely the same elements as we have derived earlier for TDHF. Completing the rest, we have:
\(R(\omega) = M(\omega)^{-1} = \begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \\ \end{bmatrix} - \hbar\omega \begin{bmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{-1} \\ \end{bmatrix} \ \ \ \ \ (34)\)

With \({\mathbf{B} = \langle ab \vert \vert ij \rangle}\). Now, we saw from the definition of the resolvent that \({R(\omega) \rightarrow \infty}\) at \({\omega \rightarrow \omega_{0n}}\), which are the poles of \({R(\omega)}\). Therefore, \({R(\omega)^{-1} \rightarrow 0}\) and
\(\begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \\ \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \\ \end{bmatrix} = \hbar\omega \begin{bmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{-1} \\ \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \\ \end{bmatrix} \ \ \ \ \ (35)\)

which are the RPA/TDHF equations. The eigencolumns determine the linear combinations of excitation and de-excitation operators that produce corresponding approximate excited states when working on the reference state \({\vert \psi_0 \rangle}\).

Derivation of general polarization propagator methods

Most higher-order response methods use some form of the polarization propagator, which is what I intend to derive here. It can then be used to derive other methods, such as TDHF/RPA and SOPPA.

We have derived earlier a response function, or frequency dependent polarizability,

where \({A}\) is the applied perturbation, and \({B}\) is the observable, and both are assumed to be Hermitian. \({\omega_{0n}}\) is the excitation energy for the change between states \({0}\) and \({n}\). It should be clear that the response function has poles when \({\omega}\) — the applied field frequency – equals to the excitation energy \({\omega_{0n}}\). Finding these poles is precisely the goal of polarization propagator methods. In the polarization propagator approach, the above equation has \({\eta}\) set to 0, and the response function (the `propagator’), defined as:
\(\langle \langle B;A \rangle \rangle_{\omega} \equiv \sum\limits_{n\neq0} \left\{ \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega - \hbar\omega_{0n}} + \frac{\langle 0 \vert A \vert n\rangle \langle n \vert B \vert 0 \rangle}{-\hbar\omega - \hbar\omega_{0n}} \right\} \ \ \ \ \ (2)\)

Now we want to describe the propagator in terms of commutators between \({A}\) and \({B}\). Make the observation that \({\frac{ab}{c+d} = \frac{ab}{c} - \frac{d}{c}\left(\frac{ab}{c+d}\right)}\), and applying to the first term of the above yields:
\(\sum\limits_{n\neq0} \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega - \hbar\omega_{0n}} = \sum\limits_{n\neq0} \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega} -\sum\limits_{n\neq0} \frac{-\hbar\omega_{0n}\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega\left(\hbar\omega - \hbar\omega_{0n}\right)} \ \ \ \ \ (3)\)

Do the same for the second term and combine, recognizing that the \({n=0}\) term vanishes in the first part (thus we get a sum over all \({n}\)), and making use of the fact that \({1 = \sum\limits_{n}\vert n\rangle\langle n\vert }\) and \({H\vert n \rangle = E_n\vert n\rangle}\) and \({\hbar\omega_{0n} = E_n - E_0}\):
\(\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \langle 0 \vert \left[B,A\right] \vert 0 \rangle + \frac{1}{\hbar\omega} \sum\limits_{n\neq0} \left\{\frac{\langle 0 \vert B \vert n\rangle\langle n \vert \left[H,A\right] \vert 0\rangle}{\hbar\omega - \hbar\omega_{0n}} + \frac{\langle 0 \vert \left[H,A\right] \vert n\rangle\langle n \vert B \vert 0\rangle}{-\hbar\omega - \hbar\omega_{0n}}\right\} \ \ \ \ \ (4)\)

Which is to say that
\(\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \langle 0 \vert \left[B,A\right] \vert 0 \rangle + \frac{1}{\hbar\omega}\langle \langle B;\left[H,A\right] \rangle \rangle_{\omega} \ \ \ \ \ (5)\)

Or, as we will use it:
\(\hbar\omega\langle \langle B;A \rangle \rangle_{\omega} = \langle 0 \vert \left[B,A\right] \vert 0 \rangle - \langle \langle \left[H,B\right] A \rangle \rangle_{\omega} \ \ \ \ \ (6)\)

As you may have started to see, we can define the propagator iteratively in terms of commutator expectation values of ever-increasing complexity. This is what is known as the so-called ‘‘moment expansion’’ of the propagator. Thus by iteration:
\(\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \left\{ \langle 0 \vert \left[B,A\right] \vert 0 \rangle + \left(\frac{-1}{\hbar\omega}\right)\langle 0 \vert \left[\left[H,B\right],A\right] \vert 0 \rangle + \left(\frac{-1}{\hbar\omega}\right)^2\langle 0 \vert \left[\left[H,\left[H,B\right]\right],A\right] \vert 0 \rangle + \cdots \right\} \ \ \ \ \ (7)\)

We introduce the ‘‘superoperator’’ (analogous to the Liouville operator in Statistical Mechanics), which acts on operators to give their commutator:

With this definition, we have the power series
\(\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \sum\limits_{n=0}^{\infty} \left(\frac{-1}{\hbar\omega}\right)^n \langle 0 \vert \left[\hat{H}^nB,A\right]\vert 0\rangle \ \ \ \ \ (9)\)

At this point we make two useful observations. First, recognize that
\(\langle 0 \vert \left[\hat{H}B,A\right]\vert 0\rangle = -\langle0\vert \left[B,\hat{H}A\right]\vert 0\rangle \ \ \ \ \ (10)\)

and so \({\hat{H}}\) can be applied to \({A}\) instead of \({B}\) insofar as we introduce a factor of \({(-1)^n}\). Furthermore, note that the power series is equivalent to
\(\frac{1}{1-x} = 1 + x + x^2 + x^3 + \cdots \ \ \ \ \ (11)\)

Making use of these two observations (and using \({\hat{1}X = X}\) and \({\hat{H}^0 = \hat{1}}\), where \({\hat{1}}\) is the unit superoperator), we have
\(\langle \langle B;A \rangle \rangle_{\omega} = \langle 0 \vert \left[B,\left(\hbar\omega\hat{1} - \hat{H}\right)^{-1}A\right]\vert 0\rangle \ \ \ \ \ (12)\)

Which is merely a cosmetic change at this point, as the superoperator resolvent is defined by the series expansion. We need to find a matrix representation of the resolvent, which implies that we find a complete basis set of operators. To do this, we are going to develop an operator space, where \({\hat{H}}\) is defined by its effect on operators instead of vectors. Introducing the notation
\(\left(X\vert Y\right) = \langle 0 \vert \left[X^{\dagger},Y\right] \vert 0 \rangle \ \ \ \ \ (13)\)

and it follows that \({\left(Y\vert X\right) = \left(X\vert Y\right)^*}\). As defined, we now have
\(\langle \langle B;A \rangle \rangle_{\omega} = \left(B^{\dagger}\vert \left(\hbar\omega\hat{1} - \hat{H}\right)^{-1}\vert A \right) \ \ \ \ \ (14)\)

Which is formally exact, albeit useless until we develop approximations. However, the form of the above equation does look similar to ordinary vector spaces in Hartree-Fock, etc. methods. Truncation of a basis in linear vector space \({V}\) to \({n}\) elements produces a subspace \({V_n}\), and truncation of a general vector corresponds to finding its projection onto the subspace. It follows, then, that we need to find a projection operator \({\rho}\), associated with the truncated basis. If the basis (\({\mathbf{e}}\), say) is orthonormal we write
\(\rho = \sum\limits_i e_ie_i^* = \mathbf{ee}^{\dagger} \ \ \ \ \ (15)\)

which in a complete basis gives:
\(\rho = \sum\limits_i \vert e_i\rangle \langle e_i\vert = \mathbf{1} \ \ \ \ \ (16)\)

If it is not an orthonormal basis, we must include the metric matrix \({\mathbf{S} = \mathbf{e}^{\dagger}\mathbf{e}}\) (or L"{o}wdin basis \({\mathbf{\overline{e}}\mathbf{S}^{-1/2}}\)):
\(\rho = \mathbf{eS}^{-1}\mathbf{e}^{\dagger} = \sum\limits_{ij} e_i (S^{-1})_{ij}e_j^* \ \ \ \ \ (17)\)

When using a truncated basis in operator space, two kinds of projections are useful (Löwdin, 1977, 1982),
\(A' = \rho A \rho, \qquad A'' = A^{1/2} \rho A^{1/2} \ \ \ \ \ (18)\)

which are the outer projection and inner projection, respectively, onto space \({V_n}\) defined by \({\rho}\). Note that \(AB = C\) does not imply \(A'B' = C'\) does not imply \(A''B'' = C''\). Plugging the metric into \({A''}\):
\(A'' = A^{1/2}\mathbf{eS}^{-1}\mathbf{e}^{\dagger}A^{1/2} \ \ \ \ \ (19)\)

and we define
\(\mathbf{f} \equiv A^{1/2}\mathbf{e} = \left(A^{1/2}\mathbf{e}_1 \quad A^{1/2}\mathbf{e}_1 \quad \cdots \right) \ \ \ \ \ (20)\)

We assume that \({A}\) is Hermitian and positive-definite, so that \({A^{1/2}}\) can be defined. Note that \({\mathbf{S} = \mathbf{e}^{\dagger}\mathbf{e} = \left(A^{-1/2}\mathbf{f}\right)^{\dagger}\left(A^{-1/2}\mathbf{f}\right) = \mathbf{f}^{\dagger}A^{-1}\mathbf{f} \implies A'' = \mathbf{f}\left(\mathbf{f}^{\dagger}A^{-1}\mathbf{f}\right)^{-1}\mathbf{f}^{\dagger}}\). Because \({A}\) is arbitrary, replace it with \({A^{-1}}\), and since \({\mathbf{f}^{\dagger}A\mathbf{f} = \mathbf{A}}\) with \({A_{ij} = \langle f_i \vert A\vert f_j\rangle}\):
\(\left(\mathbf{A}^{-1}\right)'' = \mathbf{f}\left(\mathbf{f}^{\dagger}A\mathbf{f}\right)^{-1}\mathbf{f}^{\dagger} = \mathbf{fA}^{-1}\mathbf{f}^{\dagger} \ \ \ \ \ (21)\)

As the basis \({V_n \rightarrow V}\), the inner projection \({\rightarrow \mathbf{A}^{-1}}\), else it is simply a finite basis approximation to the inverse. This is the operator inverse in terms of a matrix inverse. Since \({\mathbf{e}}\) was an arbitrary basis defining \({V_n}\), let \({\mathbf{f}}\) define n-dimensional subspace \({V_n'}\). Thus:
\(\mathbf{A}^{-1} \approx \mathbf{eA}^{-1}\mathbf{e}^{\dagger} = \sum\limits_{ij} e_i(\mathbf{e}^{\dagger}\mathbf{Ae})^{-1}_{ij}\mathbf{e}^*_j \ \ \ \ \ (22)\)

Thus the inner projection leads to an approximation for the projector. Let us define the (as of yet undefined) operator basis:
\(\mathbf{n} = (\mathbf{n}_1 \quad \mathbf{n}_2 \quad \mathbf{n}_3 \quad \cdots) \ \ \ \ \ (23)\)

Given that the binary product (compare with \({\langle x\vert y\rangle = x^*y}\) for vectors)
\(X^{\dagger} \cdot Y = (X\vert Y) = \langle 0 \vert \left[X^{\dagger},Y\right] \vert 0 \rangle \ \ \ \ \ (24)\)

then for our resolvent superoperator we have
\(\hat{R}(\omega) = (\hbar\omega\hat{1} - \hat{H})^{-1} = \mathbf{n}R(\omega)\mathbf{n}^{\dagger} = \sum\limits_{r,s} \mathbf{n}_r[R(\omega)]_{rs}\mathbf{n}_{s}^{\dagger} \ \ \ \ \ (25)\)

where \({\mathbf{n}_r}\) and \({\mathbf{n}_s^{\dagger}}\) are the analogues of \({\mathbf{e}}\) and \({\mathbf{e}^*}\) in operator space. Finally, if
\(R(\omega) = M(\omega)^{-1}, \qquad M(\omega) = \mathbf{n}^{\dagger} (\hbar\omega\hat{1} - \hat{H})\mathbf{n} \ \ \ \ \ (26)\)

then we have
\(\langle \langle B;A \rangle \rangle_{\omega} = B\cdot R(\omega)\cdot A = B \cdot \mathbf{n} R(\omega) \mathbf{n}^{\dagger} \cdot A \ \ \ \ \ (27)\)

which is the key to calculating approximations to response properties. The matrix M is determined once we have chosen an operator basis. This approximation depends on two things: 1) the basis \({\mathbf{n}}\) (and its truncations), and 2) the reference function, that is not the exact ground state. Any approximations to these two things are where we get out various response methods.

MP2 with Python

If you have the transformed two electron integrals after a Hartree-Fock calculation, it’s trivial to perform Moller-Plesset Second Order Perturbation theory (MP2). It is a non-iterative way to calculate electron correlation energy (e.g. the energy obtained from having electrons “avoid” each other in molecules). The expression for MP2 energy contribution (in a spin orbital basis) looks like this:

Where i and j are the indices of the occupied orbitals, and a and b are the indices of the virtual orbitals. The numerator is a “double-bar integral”, which accounts for the Coulomb and Exchange energies between pairs of electrons, and the bottom are the orbital energies (eigenvalues). The numerator is a result of transforming the two-electron integrals to a molecular orbital basis, and the orbital energies are the result of a Hartree-Fock calculation. (Well, strictly speaking all of this is from a Hartree-Fock calculation…) You can find a derivation of this in Szabo and Ostlund (1996).

This is actually one of the easiest electronic structure methods to program. It is non-iterative (which means no guess vectors, chaotic behavior, etc), and can be readily constructed from integrals on hand.

Computationally speaking, we can see that we must construct four loops, to loop over our four indices. The first two must loop over occupied orbitals, and the second two must loop over virtual orbitals. I’ve attached some Python code that accomplishes this. Similar to the CIS and TDHF code, I have hard coded the transformed integrals for HeH+ in an STO-3G basis, at a bond length of 0.9295 Angstroms.

#!/usr/bin/python

####################################  
#  
#  Moller-Plesset Second Order  
#   Perturbation Theory (MP2)  
#  
####################################

from __future__ import division  
import math  
import numpy as np

####################################  
#  
#   FUNCTIONS  
#  
####################################

# Return compound index given four indices  
def eint(a,b,c,d):  
    if a > b: ab = a*(a+1)/2 + b  
    else: ab = b*(b+1)/2 + a  
    if c > d: cd = c*(c+1)/2 + d  
    else: cd = d*(d+1)/2 + c  
    if ab > cd: abcd = ab*(ab+1)/2 + cd  
    else: abcd = cd*(cd+1)/2 + ab  
    return abcd

# Return Value of spatial MO two electron integral  
# Example: (12\vert 34) = tei(1,2,3,4)  
def teimo(a,b,c,d):  
    return ttmo.get(eint(a,b,c,d),0.0e0)

####################################  
#  
#  INITIALIZE ORBITAL ENERGIES  
#  AND TRANSFORMED TWO ELECTRON  
#  INTEGRALS  
#  
####################################

Nelec = 2 # we have 2 electrons in HeH+  
dim = 2 # we have two spatial basis functions in STO-3G  
E = [-1.52378656, -0.26763148]  
ttmo = {5.0: 0.94542695583037617, 12.0: 0.17535895381500544, 14.0: 0.12682234020148653, 17.0: 0.59855327701641903, 19.0: -0.056821143621433257, 20.0: 0.74715464784363106}

####################################################  
#  
#  CONVERT SPATIAL TO SPIN ORBITAL MO  
#  
####################################################

# This makes the spin basis double bar integral (physicists' notation)  
# We double the dimension (each spatial orbital is now two spin orbitals)  
dim = dim*2  
ints=np.zeros((dim,dim,dim,dim))  
for p in range(1,dim+1):  
    for q in range(1,dim+1):  
        for r in range(1,dim+1):  
            for s in range(1,dim+1):  
                value1 = teimo((p+1)//2,(r+1)//2,(q+1)//2,(s+1)//2) * (p%2 == r%2) * (q%2 == s%2)  
                value2 = teimo((p+1)//2,(s+1)//2,(q+1)//2,(r+1)//2) * (p%2 == s%2) * (q%2 == r%2)  
                ints[p-1,q-1,r-1,s-1] = value1 - value2

#####################################################  
#  
#  Spin basis fock matrix eigenvalues  
#  
#####################################################

fs = np.zeros((dim))  
for i in range(0,dim):  
     fs[i] = E[i//2]  
    #fs = np.diag(fs)

######################################################  
#  
#  MP2 Calculation  
#  
######################################################

# First two loops sum over occupied spin orbitals  
# which is equal to the number of electrons. The last  
# two loops loop over the virtual orbitals.  
EMP2 = 0.0  
for i in range(0,Nelec):  
    for j in range(0,Nelec):  
        for a in range(Nelec,dim):  
            for b in range(Nelec,dim):  
                EMP2 += 0.25*ints[i,j,a,b]*ints[i,j,a,b]/(fs[i] +fs[j] -fs[a] - fs[b])

print "E(MP2) Correlation Energy = ", EMP2, " Hartrees"

If you run the code you’ll get a correlation energy of -0.00640 Hartrees. This is pretty small, and in fact most correlation energies are. Hartree-Fock does a good job at recovering ~99% of the energy of a molecule. However, correlation energies still translate to the order of 10 kcal/mol, which is chemically significant. This is why including correlation effects are so important. It’s hard to make accurate predictions about reactivities with such large error bars when you neglect correlation. MP2 is an easy way to do so: much more accurate methods, like coupled-cluster methods, will give you energetics to chemical accuracy (chemical accuracy is ~0.1 kcal/mol). However, they are often very costly!