Derivation of Time Dependent Hartree-Fock (TDHF) Equations

This derivation is largely based off the TD-DFT derivation (TDDFT and TDHF are very similar), found in the wonderful review, “Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules” by Dreuw and Head-Gordon (Chem Reviews, 2005).

We start with the time-dependent Schrödinger equation, given in Hartree-Fock form, along with its adjoint:

Now,

Substituting the first two expressions into the right hand side, these two expressions yields the Dirac form of the time-dependent Hartree-Fock equations,

where \({\mathbf{F}}\) is the Fock matrix and \({\mathbf{P} = \mathbf{CC}^{\dagger}}\) is the density matrix, same as in the time independent Hartree-Fock equations. Before a time dependent perturbation is applied, we assume that the system is in its electronic ground state, such that

and

with \({\mathbf{P}^{(0)}}\) and \({\mathbf{F}^{(0)}}\) as the unperturbed density and Fock matrix, respectively. The first condition comes from the fact that commuting matrices share a common set of eigenvectors, and the second comes from the fact that at convergence the eigenvectors of the Fock matrix are orthonormal. (Recall that the density matrix may be written as \({\mathbf{P} = \mathbf{CC}^{\dagger}}\), where the columns of \({\mathbf{C}}\) form an orthonormal set, e.g. \({\mathbf{C}^{\dagger}\mathbf{C} = \mathbf{1}}\)). The Fock matrix elements are given by

Note the Fock matrix dependence on the density matrix. At convergence, the Fock matrix is diagonal, with the elements corresponding to the orbital energies \({\epsilon}\).

Furthermore, the density matrix enjoys the following relations:

Following convention, indices \({\{i,j,...\}}\) correspond to occupied orbitals, \({\{a,b,...\}}\) correspond to unoccupied orbitals, and \({\{p,q,...\}}\) correspond to general orbitals. In general perturbation theory, a perturbed wavefunction (or density matrix, for our purposes) can be decomposed, to first order, as

and likewise for the density-dependent Fock matrix

where the superscript (1) indicates the first-order time dependent change. Inserting the expressions 13 and 14 into the time-dependent equation 6 and collecting only the first order terms yields

The time-dependent perturbation can be described by a single Fourier component,

The matrix \({f}\) is a one-electron operator, which describes the applied perturbation. This perturbation acts on the electron density, resulting in a first order Fock matrix response

so that the overall first order change in the Fock matrix is

Similarly, the first order density response is given as

with \({d_{pq}}\) as the perturbation densities. Plugging the above expressions for \({\mathbf{F}^{(1)}}\) and \({\mathbf{P}^{(1)}}\) into the first order TDHF equation 15, and collecting the terms containing \({e^{-i\omega t}}\) gives (after some algebra):

The terms multiplied by \({e^{i \omega t}}\) give the complex conjugate of the above expression. Because the density matrix must be idempotent, we can show that

If \({\mathbf{PP} = \mathbf{P}}\), and we expand \({\mathbf{P}}\) as a perturbation with arbitrary scalar \({\lambda}\), we have:

Equating like powers of \({\lambda}\) yields the expression 21, as well as our original expression 8. Using 21 allows us to restrict the elements of \({d_{pq}}\), utilizing the nature of the zeroth order density matrix given in 11. As an example, consider the element \({d_{ii}}\):

which is only true if \({P_{ii}^{(1)} = 0}\). Similar arguments show that the only contributing blocks of \({d_{pq}}\) are \({d_{ia}}\) and \({d_{ai}}\). In other words, the only contributions to the TDHF equations are the occupied-virtual and virtual-occupied blocks. All virtual-virtual and occupied-occupied blocks are necessarily zero. With this in mind, along with the diagonal nature of the unperturbed Fock matrix, yields

and

with \({x_{ai} = d_{ai}}\) and \({y_{ai} = d_{ia}}\), as is convention. The derivatives inside the expression evaluate as follows:

Finally, we make the assumption that the perturbation is infinitely small, which sends \({f_{ai} = f_{ia} \rightarrow 0}\), and we recover (recognizing that \({F_{pp}^{(0)} = \epsilon_p}\) and \({P_{ii}^{(0)} = 1}\)):

Which is the non-Hermitian eigenvalue TDHF equation. The elements of \({\mathbf{A}}\) and \({\mathbf{B}}\) are

and

Alternate derivation of linear response function

I’ve given a derivation before of the density-density response function, but today I want to give an alternate derivation of the more general linear-response function, which will prove to be useful in the derivation of the time-dependent Hartree-Fock equations (TDHF), also known in the nuclear physics community as the Random Phase Approximation (RPA). This derivation is largely taken from McWeeny (1989).

In order to derive the response function — also called frequency-dependent polarizability — we must first partition the Hamiltonian into two parts: the time-independent Hamiltonian, and the time dependent response:

Furthermore, the time-dependent Schrodinger equation is given as:

In the interaction picture, \({\Psi'_0(t) = \Psi'_0e^{-iE_0t/\hbar}}\). Thus we can partition the time-dependent wavefunction (expanding over basis of complete eigenstates) as

Where

are real, positive, exact excitation frequencies of unperturbed system. Note that \({c_0(t) = 1}\). We are assuming that \({H'(t)}\) is turned on slowly at time \({t \rightarrow -\infty}\). Substitute 1 and 3 into 2, separate the orders, and impose the boundary conditions that \({c_0 = 1}\) and \({c_m = 0}\) \({(n\neq0)}\) at \({t \rightarrow -\infty}\). This gives

If we let

Where a ‘fixed’ Hermitian operator \({A}\) determines the `shape’ of the perturbation, while time dependence is confined to the (real) ‘strength’ factor \({F(t)}\). For a perturbation beginning at time \({t \rightarrow -\infty}\) up to time \({t}\),

Which, to first order, determines the perturbed wavefunction. Now we are interested not in the perturbed wavefunction per se, but rather in the response of an observable \({O}\) to the perturbation.

where

This is a time correlation function, relating fluctuation of \({\langle O \rangle}\) at time t to the strength of the perturbation \({A}\) at some earlier time \({t'}\). \({K(OA\vert t-t')}\) is defined only for \({t'<t}\), in accordance with the principle of causality. Thus, it is a function only of the difference \({\tau = t - t'}\). Recalling the definitions of the Fourier transform \({f(\omega)}\):

Then instead of 8, we have:

Requiring \({H'(t)}\) to be Hermitian,

Now,

Which, upon combining the expressions for \({H'(t)}\) so as to `Hermitize’ the expression:

Thus

with \({F(t)}\) real. Instead of working in the time domain, we may also consider the response in terms of a single oscillatory perturbation. This means that

To ensure \({H'(\omega)}\) builds smoothly from zero at \({t \rightarrow -\infty}\), we can introduce a convergence factor \({e^{\eta t}}\) with the initial condition \({c_0 = 1}\) and \({c_n = 0}\), which gives:

Then, collecting terms of \({\pm \omega}\):

Finally:

Which is the response function, or frequency-dependent polarizability.

Linear Response TD-DFT: The density-density response function

Linear response time-dependent density function theory (LR-TDDFT) is kind of a mouthful to say, but is a really simple idea in practice. In computational chemistry, many ideas are quite simple, but they are so laden with jargon it just scares other people away! LR-TDDFT is one way of determining single excitation energies of a molecule from first principles. It is a way of describing how matter interacts with light (or any electromagnetic field). An everyday example is color: sunlight is absorbed by some material, which excites an electron. When this electron finally relaxes, it emits light of a certain energy – that light is the color you see!

To derive the LR-TDDFT working equations (e.g. the equations you can implement on a computer), we start by deriving a density-density response function. This function describes the response of electron density (the ‘density’ of ‘density functional theory’) to an external field: light of some kind, or a magnetic field, or whatever. The potential and the electronic density ‘couple’ in some way.

We describe this interacting system with a partitioned Hamiltonian:

Where \(H_0\) is time-independent unperturbed system, and \(V_{ext}\) is the time-_dependent_ Hamiltonian describing the field potential. Investigating the potential a little more closely, we find that the potential is:

Where the density operator \(p(\textbf{r})\) is defined as:

Now, we will consider the time-dependent Schrodinger equation, but we will use the interaction picture, which describes the interacting system’s time dependence as:

With \(\vert \Psi(t)_I \rangle = e^{iH_0t}\vert \Psi(t)_I\rangle\) and \(V_I = e^{iH_0t}V_{ext}e^{-iH_0t}\).

Now, we can solve this interaction picture as a Dyson series, which we truncate after the first-order term (this is where we get linear response theory). In principle, we could keep going to get quadratic terms and so on, but the equations get much more complicated. And so:

A few notes about this expression are in order. The first term is simply the system at time zero. This is our first approximation. The second term introduces the perturbation of the potential, acting on the unperturbed wavefunction. The lower limit of infinity means that we “turn on” the perturbation infinitely long ago: this is an adiabatic approximation. In other words, we turn on the perturbation slowly enough that the ground state wavefunction can track with the disturbance. If we turn on the perturbation too sharply, the response won’t be smooth, and will display complicated ringing behavior (which we don’t care about anyway). This equation is also only valid for small perturbations – an approximation that holds true in most situations!

Now, we don’t care too much about the time evolution of the wavefunction per se. What we really care about is the response of some observable \(O\) (energy, whatever). In general, we can write this response as:

Note that the term \([O_I,V_I]\) is just the commutator expression. The \(e^{\eta t'}\) ensures adiabatic switching on, and \(\eta\) is an infinitesimal real number. Eventually we will send \(\eta\) to zero, but for now we need it for \(t \to 0\) compared with \(-\infty\). If you still don’t buy that explanation, consider what the function does at \(t = 0\) and \(t = -\infty\): at time \(t = 0\) the exponential is one, at time \(t = -\infty\), the exponential is zero.

With that said, we need to derive an expression for \(V_I\). Since \(V_{ext}\) commutes with \(H_0\), we have:

Where \(\rho_I(\mathbf{r}) = e^{iH_0t}\rho(\mathbf{r})e^{-iH_0t}\). If we substitute this expression into our expression for \(\delta\langle O\rangle\), and perform the integration from \(t' = -\infty\) to \(t' = t\), we get (after a little math):

Where \(\chi(\mathbf{r},t;\mathbf{r}',t')\) is our response function, defined as:

The term \(\theta(t - t')\) is the Heaviside function. This is a step function which is zero for \(t - t' < 0\) and 1 for \(t - t' > 0\). In other words, there is no response before we turn on the perturbation, and a full response once we do. We keep this to preserve causality of the response (no sense responding before there is something to respond to!)

Quick recap of what we have done so far. At this point we have a function for the response of any observable to some external potential (a laser, a magnet, whatever). In the case of LR- TD-DFT, we are concerned with the density response to an external potential. And what do you know? Density is an observable! In fact, we gave the operator earlier in our derivation. So the next step is to plug in the density operator for \(\langle O\rangle\) and carry the results through. Plugging in we have:

with

\(\chi (\mathbf{r},t;\mathbf{r}',t') = \lim_{\eta \to 0} -i\theta(t-t')\langle\Psi(0)\vert [\rho_I(\textbf{r},t),\rho_I(\textbf{r}',t')]\vert \Psi(0)\rangle e^{\eta(t'-t)}\).

Nothing too fancy here, just a bit of math. If we use our definition of \(\rho_I(\textbf{r},t) = e^{iH_0t}\rho(\textbf{r})e^{-iH_0t}\), the response function can be written as:

Here is where it gets interesting. Up until now, our response function has been in the time domain. If we introduce a Fourier transform, and make use of the convolution theorem, we can transform the time-dependent equations to the equivalent frequency-dependent equations. This allows us to determine responses to individual frequencies of light. This is perfect, because it allows for us to later determine specific energies (i.e frequencies) of light that give distinct electronic excitations. Performing the aforementioned Fourier transform and convoluting, we get:

with

Alright, let’s wrap this up and get the spectral representation of the density response function. If the system is initially in the ground state, e.g. ( \(\vert \Psi(0)\rangle = \vert \Psi_0\rangle\) with lowest eigenvalue \(E_0\) ), and \(H_0\vert \Psi_n\rangle = E_n\vert \Psi_n\rangle\), we can insert the complete set of unperturbed states ( \(\sum\limits_n \vert \Psi_n\rangle\langle\Psi_n\vert\) = 1 ) into the response function to get:

And there you have it! The spectral representation of the response function. The function has poles when \(\omega \to (E_n - E_0)\), which are the excitation energies. Thus we have a way to determine the excitation energies of the system! Later, I will show how we can use this result in the Kohn-Sham scheme to get excitation energies from the density functional theory framework.

Hartree-Fock Self Consistent Field Procedure

Quick links to the program and files you’ll need

SCF program: scf.py

Overlap integrals: s.dat,

Kinetic Energy integrals:t.dat,

Potential Energy integrals: v.dat,

Nuclear Repulsion: enuc.dat,

Electron Repulsion integrals: eri.dat

I do quite a bit of development work as a graduate student, and one of the most helpful parts of my education was writing a Hartree-Fock SCF (self-consistent field) program from scratch. The SCF procedure is the workhorse of most computational chemistry software packages, as it is generally the first step before doing more advanced calculations (such as MP2, etc). It is also the conceptual basis for molecular orbital theory. I’ve shared my code below, in case you want to try it out for yourself. I wrote it in Python 2.7 (it won’t work in Python 3), utilizing the NumPy package. I find it to be very readable, which makes testing out ideas fast and (relatively) intuitive. The program is written to give an idea how SCF calculations work, so it isn’t that efficient – but I only use it for two electron systems anyway so who cares!

#!/usr/bin/python

####################################  
#  
# SELF CONSISTENT FIELD METHOD  
#  
####################################

from __future__ import division  
import sys  
import math  
import numpy as np  
from numpy import genfromtxt  
import csv

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

# Symmetrize a matrix given a triangular one  
def symmetrize(a):  
    return a + a.T - np.diag(a.diagonal())

# Return compund 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 two electron integral  
# Example: (12\vert 34) = tei(1,2,3,4)  
# I use chemists notation for the SCF procedure.  
def tei(a,b,c,d):  
    return twoe.get(eint(a,b,c,d),0.0)

# Put Fock matrix in Orthonormal AO basis  
def fprime(X,F):  
    return np.dot(np.transpose(X),np.dot(F,X))

# Diagonalize a matrix. Return Eigenvalues  
# and non orthogonal Eigenvectors in separate 2D arrays.  
def diagonalize(M):  
    e,Cprime = np.linalg.eigh(M)  
    #e=np.diag(e)  
    C = np.dot(S_minhalf,Cprime)  
    return e,C

# Make Density Matrix  
# and store old one to test for convergence  
def makedensity(C,P,dim,Nelec):  
    OLDP = np.zeros((dim,dim))  
    for mu in range(0,dim):  
        for nu in range(0,dim):  
            OLDP[mu,nu] = P[mu,nu]  
            P[mu,nu] = 0.0e0  
            for m in range(0,Nelec//2):  
                P[mu,nu] = P[mu,nu] + 2*C[mu,m]*C[nu,m]  
    return P, OLDP  
# Make Fock Matrix  
def makefock(Hcore,P,dim):  
    F = np.zeros((dim,dim))  
    for i in range(0,dim):  
        for j in range(0,dim):  
            F[i,j] = Hcore[i,j]  
                for k in range(0,dim):  
                    for l in range(0,dim):  
                        F[i,j] = F[i,j] + P[k,l]*(tei(i+1,j+1,k+1,l+1)-0.5e0*tei(i+1,k+1,j+1,l+1))  
    return F

# Calculate change in density matrix  
def deltap(P,OLDP):  
    DELTA = 0.0e0  
    for i in range(0,dim):  
        for j in range(0,dim):  
            DELTA = DELTA+((P[i,j]-OLDP[i,j])**2)  
    DELTA = (DELTA/4)**(0.5)  
    return DELTA

# Calculate energy at iteration  
def currentenergy(P,Hcore,F,dim):  
    EN = 0.0e0  
    for mu in range(0,dim):  
        for nu in range(0,dim):  
            EN = EN + 0.5*P[mu,nu]*(Hcore[mu,nu] + F[mu,nu])  
    return EN

####################################  
#  
# FORM CORE HAMILTONIAN  
#  
####################################  
Nelec = 2 # The number of electrons in our system

# ENUC = nuclear repulsion, Sraw is overlap matrix, Traw is kinetic energy matrix,  
# Vraw is potential energy matrix

ENUC = genfromtxt('./enuc.dat',dtype=float,delimiter=',')  
Sraw = genfromtxt('./s.dat',dtype=None)  
Traw = genfromtxt('./t.dat',dtype=None)  
Vraw = genfromtxt('./v.dat',dtype=None)

# dim is the number of basis functions  
dim = int((np.sqrt(8*len(Sraw)+1)-1)/2)

# Initialize integrals, and put them in convenient Numpy array format  
S = np.zeros((dim,dim))  
T = np.zeros((dim,dim))  
V = np.zeros((dim,dim))

for i in Sraw: S[i[0]-1,i[1]-1] = i[2]  
for i in Traw: T[i[0]-1,i[1]-1] = i[2]  
for i in Vraw: V[i[0]-1,i[1]-1] = i[2]

# The matrices are stored triangularly. For convenience, we fill  
# the whole matrix. The function is defined above.  
S = symmetrize(S)  
V = symmetrize(V)  
T = symmetrize(T)

Hcore = T + V  
#print Hcore

#####################################  
#  
# TWO ELECTRON INTEGRALS  
#  
#####################################

# Like the core hamiltonian, we need to grab the integrals from the  
# separate file, and put into ERIraw (ERI = electron repulsion integrals).  
# I chose to store the two electron integrals in a python dictionary.  
# The function 'eint' generates a unique compund index for the unique two  
# electron integral, and maps this index to the corresponding integral value.  
# 'twoe' is the name of the dictionary containing these.

ERIraw = genfromtxt('./eri.dat',dtype=None)  
twoe = {eint(row[0],row[1],row[2],row[3]) : row[4] for row in ERIraw}

#####################################  
#  
# SCF PROCEDURE  
#  
#####################################

# Step 1: orthogonalize the basis (I used symmetric orthogonalization,  
# which uses the S^(-1/2) as the transformation matrix. See Szabo and Ostlund  
# p 143 for more details.

SVAL, SVEC = np.linalg.eig(S)  
SVAL_minhalf = (np.diag(SVAL**(-0.5)))  
S_minhalf = np.dot(SVEC,np.dot(SVAL_minhalf,np.transpose(SVEC)))  
# This is the main loop. See Szabo and Ostlund p 146.  
# Try uncommenting the print lines to see step by step  
# output. Look to see the functions defined above for a better  
# understanding of the process.

P = np.zeros((dim,dim)) # P is density matrix, set intially to zero.  
DELTA = 1.0  
convergence = 0.00000001  
G = np.zeros((dim,dim)) # The G matrix is used to make the Fock matrix  
while DELTA > convergence:  
    F = makefock(Hcore,P,dim)  
   # print "F = \n", F  
    Fprime = fprime(S_minhalf,F)  
   # print "Fprime = \n", Fprime

    E,Cprime = np.linalg.eigh(Fprime)  
    C = np.dot(S_minhalf,Cprime)  
    E,C = diagonalize(Fprime)  
   # print "C = \n", C

    P,OLDP = makedensity(C,P,dim,Nelec)

   # print "P = \n", P

# test for convergence. if meets criteria, exit loop and calculate properties of interest

    DELTA = deltap(P,OLDP)  
   #print "E= ",currentenergy(P,Hcore,F,dim)+ENUC  
   #print "Delta = ", DELTA,"\n"

    EN = currentenergy(P,Hcore,F,dim)  
    print "TOTAL E(SCF) = \n", EN + ENUC  
   #print "C = \n", C

To use it, download the scf.py file to your computer, and put it in a folder along with s.dat, t.dat, v.dat, and enuc.dat. These are the integral values it needs to run: s.dat is the overlap matrix values, t.dat is the kinetic energy matrix, v.dat is the potential energy matrix, and enuc.dat contains the nuclear repulsion energy. You’ll also need the eri.dat, which contains the two electron integrals. All values were taken from Gaussian09 for HeH+ at a bond length of 0.9295 Angstrom with an STO-3G basis set. Once you have those downloaded, fire up your terminal and give it a try! See the example of my terminal below for usage. Just navigate to the folder you have everything in, and type python scf.py then hit Enter!