Hartree—Fock Approximate Molecular Orbital Theory Justin T. Fermann 3 lectures on theory 1 lecture on programming Purpose First, we define the problem, beginning with the Schrodinger equation = EV. (1) Our goal is to come up with an analytic equation for the energy which can be minimized with respect to some variational parameter(s) to give an upper bound on the energy. To do this, we must 1. Understand the hamiltonian operater H. 2. Find an appropriate wave function \I> which allows simple calculation of the electronic energy. 3. Examine potential variational parameters to figure out how to minimize the energy. 1 It will likely be convenient to have an outline of where we're going, so here it is: I. Puropse (which we've been through) II. The Hamiltoniaii, H • Kinetic energy operators • Coulombic potential operators III. The Wave Function, * • 1 electron orbitals • Hartree products • Slater determinants IV. Hamiltonian as Energy Operator, Schrödinger Equation • E = (V\H\V) • Integrals over one electron operators • Integrals over two electron operators • Specific case energy expressions • General form of HF energy V. What Variational Parameter? • LCAO-MO theory, energy in AO basis • Density Matrices VI. Hartree-Fock Equations • Lagrangs's Undetermined Multipliers • A load 'o math VII. Matrix Formalism VIII. Program Outline 2 The Hamiltonian, H The Hamiltonian is the total energy operator for a system, and is written as the sum of the kinetic energy of all the components of the system and the internal potential energy. In an atom or molecule, comprised of positive nuclei and negative electrons, the potential energy is simply that due to the coulombic interactions present. Thus for the kinetic energy in a system of M nuclei and N electrons: M ■ N N , And for the potential energy: Te = -E^- (3) vNN = J2- (4) A>B N ! TAB Vee = E— (5) i>j r M N y Ve* = -EE^- (6) A i 'iA Since H = T + V, i 1 A i iyj rl3 A IMA A>B TAB Within the Born-Oppenheimer approximation, we assume the nuclei are held fixed while the electrons move really fast around them, (note: Mp/Me « 1840.) In this case, nuclear motion and electronic motion are seperated. The last two terms can be removed from the total hamiltonian to give the electronic hamiltonian, He, since Vjvjv = K, and = 0. The nuclear motion is handled in a rotational/vibrational analysis. We will be working within the B-0 approximation, so realizing N i M N 7 N -i 4 = -£iv?-£E^ + E^- (8) i A A i ' *A iyj 113 we completly define the problem. Solving the electronic Schrodinger equation using this will give the electronic structure of a molecular system at a fixed nuclear geometry. 3 The Wave Function, ^ We've derived a complete many-electron Hamiltonian operator. Of course, the Schrodinger equation involving it is intractable, so let's consider a simpler problem, involving the one-electron hamiltonian which involves no electron-electron interaction. This is soluable in the B-0 approximation (recall the hydrogen atom by letting M = 1). Call the solutions to the one-electron Schrodinger equation Xi(x). These will be molecular spin orbitals when we get around to it, but for now let it suffice to know they satisify the eigenequation ki)Xj(xi) = ejXjte) (10) with the interpretation that electron i occupies spin orbital Xj with energy Sj. If we ignore electron-electron interaction in He, we construct a simpler system with Hamiltonian N H = ^h(i). (11) i It will have eigenfunctions which are simple products of occupied spin orbitals, and thus an energy which is a sum of individual orbital energies, as VHP = Xi(xi)Xi(x2)x*(x3)---X,(xiv) (12) = ei + ej+ £* + ... + £„ = E. (13) This kind of wavefunction is called a Hartree Product, and it is not physically realistic. In the first place, it is an independent-electron model, and we know electrons repel each other. Secondly, it does not satisfy the antisymmetry principle due to Pauli which states that the sign of the wavefunction must be inverted under the operation of switching the coordinates of any two electrons, or *(---xi---xi---) = -*(---xi---xi---). (14) Part of the proof of equation 13 acknowledges this is not so for a Hartree Product. To remedy this, first consider a two-electron system, such as helium. Two equivalent Hartree Product wavefunctions for this system are *fp = Xi(xi)Xi(x2) ^P = X*(x2)Xi(xi). (15) 4 Obviously, neither of these is appropriate. However, using the old "by inspection..." trick, we notice that * = ^IXiMXjfa) ~ Xi(x2)Xj(x1)] (16) does. The mathematical form of this wavefunction can be generated by a determinant of x's, Xi(xi) XjM Xi(x2) Xj(x2) (17) The familiar Pauli exclusion principle follows directly from this example. When we attempt to doubly occupy a spin orbital Xi by putting electron 1 and electron 2 in it, what happens? 2-i/2 Xi(xi) Xi(xi) Xi(x2) Xi(x2) ^[Xi(xi)Xi(x2) - Xi(x2)Xi(xi)] 0 (18) Equation 17 can be generalized to give the N electron Slater determinant * = (AT!)-V2 Xi(xi) Xj(xi) ••• Xu(xi) Xi(x2) Xj(x2) ••• Xu(x2) Xi(xjv) Xj(xjv) ••• Xt)(xjv) (19) A shorthand notation for a Slater determinant has been introduced, where all the diagonal elements in the determinant are written in order as a "ket" vector. Equation 19 can thus be written as |Xi(xi)Xj(x2)---Xu(xjv)> where the normalization constant is absorbed into the notation. (20) Now we have written down a wave function appropriate for use in the case where H = Yl h{i). In HF theory, we make some simplifications so many-electron atoms and molecules can be treated this way. By tacitly assuming that each electron moves in a percieved electric field generated by the stationary nuclei and the average spatial distribution of all the other electrons, it essentially becomes an independant-electron problem. The HF Self Consistent Field procedure (SCF) will be bent on constructing each x(x) to give the lowest energy. 5 Energy Expressions Let's assume a wave function of the Slater determinant form and find an expression for the expectation value of the energy. We've written a Slater determinant as a ket vector in shorthand notation, allowing us to make use of Dirac notation for such things as overlap. In this context, recall that <*a|*6) = |*:(x)*6(x)dx (21) where the basis vectors \I> is expanded in are every possible value of x with contraction coefficients identified as the value of ^(x) at x. Thus placing an operator (such as H) inside the bracket, we get the expectation value of the observable associated with that operator. Since H is the energy operator, E0 = J dr%HeV0 = (*0\He\*0). (22) dr is the differential of all the spin and space coordinates of all the electrons. With much foresight, we continue to simplify the problem by writing H as a sum of one- and two-electron operators He = E^ + Er (23) i i>j 'U = H^ore + H2. (24) This will allow us to more precisely develop the electronic energy by it's components. First, examine the core hamiltonian H^°re. (y\Hrem = e<*i^(oi*> (25) i = E(Xi(xi) .. .Xi(xi) ■ ■ ■ IMi)IXi(xi) .. .Xi(xi) ...) (26) i The nature of this is best evidenced by example, so we turn to the familiar helium atom, \I> = |Xi(x1)x2(x2))- Look at one term in the above sum, for the sake of illustration take ^(1). = 1 I[x*(Xl)x;(x2) - Xt(x2)x5(x1)]^(l) [Xi(xi)x2(x2) - xi(x2)x2(xi)]dxidx2 (27) 6 = \l ^*(Xl)^2(x2)^(l)Xi(xi)x2(x2)o?r + J xl(x2)x2(xi)fc(l)xi(x2)x2(xi)dr - /xI(x2)x2(xi)A(l)xi(x1)x2(x2)dr - JxI(x1)x5(x2)A(l)xi(x2)x2(xi)dr (28) = |«Xi|^|Xi> + = (Xi|%i> + (X2|%2> (30) i Profound, isn't it? Seems that every occupied spin orbital Xi yields a term of the form (xil^lXi) to the one electron energy. Now look at H2. i>j 'V = E(Xi(xi) ■ ■ -Xi(xi) ■ ■ ■ | — |xi(xi) ■ ■ -Xi(xi) ■ ■ ■) (31) i>j riJ Continuing to work in the helium atom example (realize that this could be any two electron system) pick = (1,2) and look at that one term. 1 If 1 (Xi(xi)x2(x2)| — |xi(xi)x2(x2)) = - / [xi(xi)x2(x2) - Xi(x2)x2(xi)] — [Xi(xi)x2(x2) - Xi(x2)x2(xi)]dx!dx2 (32) If 1 f 1 = ö / *t(xi)x2(x2)—Xi(xi)X2(x2)dr + / x*(x2)x2(xi)—Xi(x2)x2(xi)dr - I xt(x2)X2(x0—Xi(xi)x2(x2)dr - / x*(xi)x2(x2)—Xi(x2)x2(xi)dr (33) •> f\2 J f\2 Unfortunately, the ^ operator prevents seperation of the integrations over the electronic coordinates of electron 1 and electron 2. It cannot be assured that the last two terms are zero. In 7 general, they are not. However, since the xi and x2 are dummy variables, the first and second terms of equation 33 are equal, as are the last two. Thus for the two electron operator (Xi(xi)x2(x2)| — |xi(xi)x2(x2)> = (12|12> - (12|21> = (12||12> (34) T12 where (ij\kl) = I dxidx2X,*(xi)Xi(x2)—X*(xi)xi(x2) (35) J n2 (ij\\kl) = (ij\kl)-(ij\lk). (36) The construction (ij\\kl) is called an antisymmetrized two electron integral in physicists notation. By working in the spin orbital basis, much trouble is avoided. In fact, by extending the results shown previously to the general case, we can now write down the HF energy for a given set of occupied spin orbitals. EHF = = (*|#iCore + H2\*) (37) Ehf = (Xi(xi) ■ ■ -Xi(xi) ■ ■ ■ |#i|xi(xi).. -Xi(xi) ■ ■ ■) + (Xi(xi) ■ ■ ■ Xi(xi) ■ ■ ■ |#2|Xi(xi) ■ ■ ■ Xi(xj) ■ ■ ■) (38) = £(Xi(xi) ■ ■ -Xi(xi) ■ ■ ■ |A(i)|xi(xi) ■ ■ -Xi(xi) ■ ■ ■) + i £(Xi(xi) ■ ■ -Xi(xj) ■ ■ ■ | —|xi(xi) ■ ■ -Xi(xj) ■ ■ ■) (39) i>j r12 = £ + £ (40) i i>j Now move on and consider working in the spatial orbital basis, where Xi(x) = ^a(r)o;. (41) This is more natural, since our intuition is usually based on having a region of space which describes the location (more or less) of two electrons, one of alpha spin and one of beta spin. Some of quantum chemistry is formulated entirely in terms of spin orbitals, for various reasons. For our purposes, we will work entirely in the spatial orbital basis. This will cause things to get somewhat murky soon, but in the long run it will be simpler. At any rate, in the two electron system we adore so much, we can identify the two occupied spin orbitals xi and X2 as the spin up and spin down halves of the single lowest lying spatial orbital, a Is in helium or the a bonding orbital in H2 for example. These can be more precisely defined as Xi(x) = ^i(r)a X2(x) = ViW. (42) (43) 8 This changes the way we write slater determinants. Using an overbar to denote (5 spin occupation of a spatial orbital, ip, |Xi(xi)x2(x2)> = \M^i(T2)) = |Vv0i> (44) Rethinking the one electron integrals for this case, Y,(Mi\kQ\Mi) = (^WI^ + ^ilMWi) (45) i = (Vi \h\ fa) + (V>iWi) (46) = /in +/in (47) The notation {t^^h^i) is used to denote an integral over only spatial coordinates, what remains after the spin integrations have been carried out, giving a factor of 1 or 0. That was a neat closed shell system. How about something like \I> = IV'iV'iV^)? £<*|fc(0l*> = (^i|^i> + (^i|^i> + (^|^> (48) i = (^i|%i) + (^i|%i) + (^|%2) (49) = 2/in + /i22 docc socc i i occ = (50) i The coefficient /j here is related to the occupation of spatial orbital i, and will be more precisely defined later. The two electron integrals are a tad more involved, but we go about it in essentially the same manner. (X1X2IX1X2) = (V'lV'llV'lV'l) (51) = WM^M (52) = ftMM^i) (53) (X1X2IX2X1) = (A^il^M (54) = [^il^iV'i] (55) = 0 (56) Yes, I know. Very confusing. But it's all just notation, and can be understood. In physicist's notation (equivalent to Dirac notation), (ipiipjlipkipi) refers to the two electron integral where tpi and tpk are functions of electron 1, while tpj and tpi are functions of electron 2. Chemist's notation (with the square brackets []) places the functions of electron 1 on the left and the functions of 9 electron 2 on the right. When the two functions of a single electron are not of the same spin, the whole integral goes to zero, otherwise the spin integrates out to 1. Hence the spin-free notation (V'iV'iIV'iV'i)- What occurs when the two electrons are of parallel spin, requiring distinct spatial orbitals and a wavefunction something like \I> = IV'iV^)? The same general form is present, and a related antisymmetrized two electron integral is evaluated. In this case, (V'iV^HV'iV^)- WiikUM = (12ll12> = (12|12) - (12|21) = [11122] - [12|21] = (11|22) - (12|21) (57) = Ji2~K12 (58) Another bit of notation, which should be apparent. Ji2 = (H|22) and K12 = (12121). Jij is termed a coulomb integral and has the physically reassuring interpretation of somehow accounting for electronic repulsion between electrons in molecular orbital i and molecular orbital j. K^, the exchange integral, has no classical analog and no true physical interpretation. Many have tried to come up with something, and a typical attempt says that it "correlates the motions of electron i and j when they have parallel spins, lowering the energy since those electrons avoid each other better." Whatever. Best to just move on to a general energy equation in the spatial MO basis. In summary: • The one electron integrals contribute ha for each electron in orbital i. • The two electron integrals contribute for each pair of electrons, and (—Kij) for each pair of parallel spin. occ occ Ehf = 2 UKi + Y^{aij(H\jj) + Pa(ij\ij)} i i,j occ occ = 2 UKi + Y,iaijJij + PijKij} (59) i i,3 fi = < 1 if i doubly occupied 1/2 if i singly occupied 2 if i and j doubly occupied 1 if i or j doubly occupied and the other singly occupied 1/2 if i and j singly occupied (60) (61) 10 -1 if i and j doubly occupied -1/2 if i or j singly occupied, the other doubly occupied -1/2 if i and j singly occupied, with parallel spin 1/2 if i and j singly occupied, with opposite spin 11 What Variational Parameter? Ah, the crux of the problem, is it not? Up until now, we've just assumed we have some set of molecular orbitals Xi or V^j which we can manipulate at will. But how does one come up with even approximate solutions to the many body Schrodinger equation without having to solve it? Start with the celebrated linear combination of atomic orbitals to get molecular orbitals (LCAO MO) approximation. This allows us to use some set of (approximate) atomic orbitals, the basis functions which we know and love, to expand the MOs in. In the most general terms, ^ = e^- (63) tpi remains a spatial molecular orbital, <^ is a spatial atomic orbital (perhaps symmetry orbital, but no matter), and C% are the contraction coefficients by which we transform from one basis to another. Armed with only this, we should be able to compose the electronic energy in the atomic orbital basis. Why, you ask? Because we have an expression in terms of integrals over MOs. To variationally minimize that energy, we need to vary the MOs themselves, but have no way to do that, since they remain these amorphous constructs. By defining them a bit more precisely, we should arrive at a point where an obvious set of variational parameters (hint: C^) present themselves. Begin with the closed shell HF energy in terms of spatial MOs. docc docc EHF = ha + ^2{2Jij — Kij} i i,j docc docc = 2Y,K + Y.Wi\n)-mJ)} Since {^i\h\^i) = J(ri)/i(l)^(ri)dri and ^ = E^CjA, /^(rOMlMCri)*! = /EC?^(n)Ml)E^(ri)*i = EEW/^(ri)Mi)^('i)*i fl V J = EEWV (65) (<^t|/i|<^) is a one electron integral over atomic orbitals. Do we have something we can actually calculate!? Take an aside and examine this quantity superficially. 12 Typically, basis functions are constructed to mimic true atomic orbitals. The Hydrogen atom can be described rigorously, and the eigenfunctions found. The Is orbital looks something like Ne~l>r. It satisfies all the appropriate boundary conditions, having a cusp at the nucleus and exponentially decaying to zero at infinity. Higher angular momentum functions, like 2p's and 3d's, can be built from this basic framework through adding the angular nodes by multiplying in factors of x, y, and z. Basis functions such as these are called slater-type orbitals. If instead of the exponential Ne~l>r we use 2 a gaussian function, Ne~ar , we loose the boundary conditions but generate a more tractable problem when it comes to calculating integrals. Using a linear combination of single cartesian gaussian-type orbitals to approximate a slater-type orbital gives better computational accuracy without too much more effort. Here's the functional forms of all three types: (j)lTO(r) = Nxlymzne-Cr (66) <^GTO(r) = Nxlymzne-ar2 (67) (71) pv pvpcr = Yl D^[2V + E D„42(HH - G"pM}] (72) p,v pa is the density matrix, a product of AO-MO coefficient matrices, or p,v = CtC (73) 13 Hartree Fock Equations The electronic energy is a functional of the spin orbitals, and we want to minimize it subject to some set of constraints. This can be done using the calculus of variations applied to functionals. So lets look at a general example of functional variation applied to E, a functional of some trial wavefunction $ that can be linearly varied under a single constraint. E = ($|#|$) (74) n \*) = E*l**> (75) i By equation 74, we see that E[$], depending on the form of the wavefunction, and by equation 75 that |$) can be linearly expanded (hence linearly varied) in some set of N functions. This is directly analogous to expanding the asymmetric top rotational wavefunctions in a complete set of symmetric top rotational wavefunctions. The task is to minimize E subject to the single constraint that the wavefunction $ remain normalized, or (m = j:^cjm^j) = i (76) ij E={*\H\*) = Y, + ^<*«l*i»] = E^[E9i«*«l^l*i> + ^l*i»] + « j E^Eci((*il^l*«> + ^<*il*«»] = E <^*E + + E o = \X1X2 • • 'Xn)- We state the problem as: Please minimize the electronic energy of this single determinant subject to the constraint that the spin orbitals all remain orthonormal to one another. 15 We already understand the energy fine by equation 37, and the constraint can be simply stated as (a\b) -Sab = 0 (85) There are N spin orbitals, so there are N(N + l)/2 independent constraints (note: (a\b) = (b\a)*), so we need that many undetermined multipliers in our lagrangian function, for which we present N a £[{Xa}} = E0[{Xa}} -ee£u(#) - S„t) a b (86) The restricted sum will prove inconvenient, but it can be eliminated. By taking the complex conjugate of the constraints and the lagrangian function, equation 86, (a\b)*-5ab = 0 (b\a)-Sab = 0 N a £[{Xa}} = ^[{X«}]-EE4,«&|fl>-$*) a b (87) (88) we realize we can unrestrict the summation by restraining eab to be a hermitian matrix, such that £ab = £ba- This introduces no new undetermined multipliers into the equation, and creates a form more amenable to further derivation. Thus the Lagrangian function we will be working with is N N £[{Xa}} = E0[{Xa}} -ee^((#) - Sat)- a b The differential of this function must be set to zero as before, giving us N N SjC = SE0-Y^2£baS(a\b) = 0 a b Since we have an electronic energy in terms of spin orbitals (89) (90) N N E0 = 5>|#|a) + -$>&||a&) a N ab N E0 = J2(a\H\a) + o J2Kab\ab) ~ (ab\ba)}, (91) ab we can write the variance of the energy SE0 as N SE0 = J2((Sa\H\a) + (a\H\Sa)) (Sab\ab) + (aSb[ — (Sab\ba) — (aSb\ba) — (ab\Sba) — (ab\bSa) 1 ^ (Sab\ab) + (aSb\ab) + (ab\Sab) + (ab\aSb)+ 2 ^ ab (92) 16 It takes some thought to realize that there are only two unique two electron integrals in this list, and that it can be written N N 5E0 = ^2(5a\H\a) + ^2((5ab\ab) — (5ab\ba)) + complex conjugate (93) a ab The other variance we need is ebaS(a\b), which can be expanded as Y,£ba(8a\b) + Y,£ba(a\Sb) = Y,£ba(Sa\b)+ Y,£ab(b\Sa) ab ab ab ab = J2(eba(6a\b) + eta(6a\by). (94) So the whole thing boils down to one neat statement, N N SC = ^2(5a\H\a) + ^2(Sa\b) + ^2((5ab\\ab) — £ba(3a\b)) + complex conjugate = 0 (95) a a ab If we conveniently define a coulomb operator Jb(l) and an exchange operator Kb(l) as Jb(l) = J d^2\xb(2)\2r^ (96) Kb(l)Xa(l) = J dx2X6*(2)rr21Xa(2)x6(l), (97) we can rewrite the two electron integral (ab\\ab) = JdxlX:(l)(Jb(l) - Kb(l))Xa(l) = (a\(Jb(l)-Kb(l))\a). (98) This allows for a more compact notation to be employed in writing the variance in the lagrangian function 6C = J2[ d^6x:(l)[h(l)xa(l) + Ei(^(l) - Kb(l))xa(l) - ebaXb(l)}] + complex conjugate = 0. a J b (99) Like before, the part in brackets is forced to be zero, since 8x*JS) can be anything. Setting it equal to zero and rearranging to make it look like some sort of eigenvalue equation yields [Ml) + EW)-^(l))]Xa(l) = E^aXfc(l) (10°) b b f(l)Xa(l)=J2£baXb(l) (101) b 17 These are the glorious Hartree Fock equations derived in general in the spin orbital basis. But wait - there's a problem. These are coupled integro-differential equations, and while they are not strictly unsolvable, they're a pain. It would be nice to at least uncouple them, so let's do that. If we apply a unitary rotation to the full set of spin orbitals, generating a new set x'a = T,XbUba (102) b where U is unitary, ie. = U-1, what changes? The rotation can be written as a matrix product, if we define A as the matrix resembling the slater determinant for the system, ie. (TV!)-1/2 det(A) = |\&0)- In that case, det(A') = det(A) det(U) (103) However, since this matrix U is unitary, | det(U)|2 = 1, det(U) = e1^ and the new wavefunction differs from the old by a phase factor, affecting nothing observable. How does it affect f(l) and E£(l) = E fdx2x':(2)^x1(2) (104) a a J = E / rfx2E KtiVYu E UcaXdl) (105) a J b c = E(E Wa) / dx2X;(2)rii1Xc(2) (106) be a J = |dx2X6*(2)r1-21Xc(2) = E ^(1) (107) be E#(l)xi(l) = E/^2x:(2)rr2^(2)xl(l) (108) a a J = E / ^X2 E ULxtWr^XdW E ^aXc(l) (109) a J b c = J2(T,ULUea) [d^xt(2)r^Xd(2)'xe(l) (HO) be a J = J26bc[ dx2X6*(2)rr21Xd(2)'xc(l) = E^(!)xi(l) (HI) be J b So /(l)' = 7(1)! How about £hQ? Start by realizing that Eba are matrix elements of the fock operator. = Ee6a (H2) b 18 = eca (113) 4 = fdxlX':f(l)x'b(l) (H4) = J2U*caU^cd (115) cd = UVU (116) This last result can be written as a martix product as well, and it is seen that this is now a unitary transformation to the matrix e. We are free to choose U to be whatever we please, and if we choose it to make e diagonal, we can rewrite the Hartree Fock equations as f\Xa) = ea\Xa). (117) When this is done, the resulting spin orbitals are termed the Hartree-Fock canonical orbitals. Read section 3.3 in Szabo and Ostlund for various fun things to do with the Hartree-Fock equations. The problem still remains, though. These are integro-differential equations, which computers (and computer programmers) balk at. That is why Roothaan is a HERO! Through his results, we can transform these into a set of matrix formulated algebraic equations that computers and programmers dig. The general case is too troublesome for now, so let's limit ourselves to closed shell, RHF orbitals. To take advantage of the simplifications this can afford, we need to return to the spatial orbital basis. We've derived the spin orbital based Fock operator N /(Xl) = h(Xl) + e(Ja(Xl) - £a(Xl)). (118) a Without further ado, we'll introduce the spatial orbital based Fock operator and be done with it. N/2 /(ri) = Mr1) + X)(2Ja(r1)-^a(r1) (119) a Jb(ri) = J dv2\^b{v2)\2r^ (120) ^(ri)Va(ri) = J dr2^(r2)rr2Va(r2)^(l) (121) With the properties /('i)^-('i) = ^i('i) (122) N/2 Si = hu + ^22Jib- Kib (123) b 19 Now we introduce a basis set expansion to bring the HF integro-differential equations to soluable algebraic equations. Letting = C^cj)^, /(i)E^ = e* (1) and integrating over electron 1 gives E<3 / ^;(i)/(i)i(i) = £iEc; / ^(1)^(1). (125) Identifying the integrals as matrix elements of the fock operator and the unit operator (overlap) respectively, E^ci^E^cj (126) V V Using the fact that e is diagonal, this can be written as the matrix product FC = SCe (127) So what is F, this so called Fock Matrix? We've defined before as the matrix element of the one electron fock operator, f(l) in equation 119. Writing out this integral and expanding, F,v = 1^^(1)7(1)^(1) = /dn^(i)[/i(i) + E(2^«(i) - J a = J ^^(1)^(1)^(1) + E[2 / dr^l) Ja(l)^(l) - J dn^(l)^a(l)^(l)] (128) = + E 2(/w|aa) — (fj,a\au) a = V + EEWH/»)-(wHl a pa = v + E(EW)[2(^lp^)-(^IH] pa a = v + ED^[2(HH-(/^M] (129) pa This is a quantity which can be easily constructed given a set of molecular orbitals (the coefficients C^) and a precalculated set of atomic orbital integrals. At this point, the Hartree Fock equations have been reduced to a matrix eigenvector problem, FC = SCe, but not in a computationally convenient form. Following the analysis leading to equation 84, we first define the transformed Fock matrix as F* = S"1/2FS-1/2. (130) 20 We can then take F^S^C) = (S^2C)e (131) to be equivalent to equation 127. Since e is diagonal by choice, diagonalizing the transformed Fock matrix gives a set of transformed molecular orbital coefficients, = S1/,2C. The matrix S1/,2C can be back transformed to give the true MO coefficient matrix, C. The density matrix for these coefficients is formed by the product D = C*C, and can subsequently be used to construct a new fock matrix via equation 129. Since the overlap matrix S does not depend on the MO coefficients, the same unitary transformation can be applied to the new fock matrix to give a new transformed fock matrix. This can be diagonalized to produce new MO coefficients, and the process repeated until convergance. As an initial guess for the fock matrix, one generally uses the core hamiltonian, ignoring all the two electron integrals. From the core hamiltonian, an initial C is obtained and a much improved Fock matrix can be built including the two electron integrals. And that's about it for the Hartree Fock Self Consistent Field Method. These last few pages will be the most important when you get around to programming the closed shell SCF method for the specific case of water, as you will be given the integrals in a file, and you can begin the process by building the core hamiltonian as described above. (132) 21