Chem 3502/5502 Physical Chemistry I (Quantum Mechanics) 3 Credits Fall Semester 2009 Laura Gagliardi Lecture 29, November 20, 2009 Solved Homework Water, H 2 O, involves 2 hydrogen atoms and an oxygen atom. To minimaly represent hydrogen is simple—one needs only a 1s function. As there are two H atoms, we wil need two such functions, one centered on each atom. Oxygen has electrons in the second principal quantum level, so we wil need one 1s, one 2s, and thre 2p functions (one each of p x , p y , and p z ). So, that's five functions total on O, plus one each on the H atoms which is then 7 total. The total number of one-electron integrals we need to solve is determined by considering that we need kinetic-energy integrals T µ! = µ " 1 2 # 2 ! and potential energy integrals ! V µ" =µ# Z k r kk nuclei $" Since both µ and ν can be any one of the 7 basis functions, there are 7 x 7 posible kinetic energy integrals (49) and the same number of potential energy integrals for each nucleus. As there are 3 nuclei, that is 7 x 7 x 3 = 147 integrals. The grand total of one- electron integrals is thus 196. The two-electron integrals are µ!"# ( ) = $ µ 1( )$ ! 1( ) 1 r 12 $ " 2( )$ # 2( )dr 1 dr 2 % where now µ, ν, λ, and σ can be any one of the seven basis functions (note that leting them be any function covers both al Coulomb and al exchange integrals). There are thus 7 x 7 x 7 x 7 = 2401 two-electron integrals. The numbers computed above involve the contracted basis functions, each of which, since the basis is STO-3G, is composed of 3 primitive functions. Thus, for any individual one-electron integral, there wil be 3 x 3 = 9 separate integrals involving the primitives. There are thus 9 x 196 = 1764 individual primitive one-electron integrals. 29-2 As for the two-electron integrals, again, every individual integral wil require considering every possible combination of constituent primitives which is 3 x 3 x 3 x 3 = 81. Thus, the total number of primitive two-electron integrals is 81 x 2401 = 194,481 (gulp!) Notice that even for this smal molecule the number of two-electron integrals totaly dominates the number of one-electron integrals. The disparity only increases with molecular size. Al this to find 5 occupied molecular orbitals from which to form a final Slater determinant (10 electrons, two to an orbital, so 5 orbitals). The situation sounds horrible, but it should be recognized that by using gausians, the solutions to all of the integrals are known to be analytic formulae involving only interatomic distances, cartesian exponents, and α values in the gausians. So, the total number of floating-point operations to solve the almost 200,000 grand-total integrals may be about 1,000,000. In computer speak that's one megaflop (megaflop = milion FLoating-point OPerations). A modern digital computer procesor can achieve gigaflop per second performance, so the computer can acomplish al these calculations in under one second. In fact, modern computers are so fast at calculations that it can be faster to recompute the integral values than to take the time to write them onto a storage device to retrieve the number later! The second way in which things can be improved is to recognize that there are symmetries of which advantage can be taken. From the turnover rule, it must be true that T µ! = µ " 1 2 # 2 ! = !" 1 2 # 2 µ =T !µ and the same is true for the nuclear atraction integrals. This reduces the total number of contracted integrals to N(N+1)/2 where N is the number of contracted basis functions. For water this wil be 28, then, instead of 49. When basis functions µ and ν are the same (but possibly on diferent atoms) there wil be additional savings from symmetry since they wil be composed of the same primitives. A similar analysis for the two-electron integrals leads to a final number of N(N+1)(N 2 +N+2)/8. For water, with N = 7, that is 406 contracted two-electron integrals—a substantial improvement over 2401. Writing computer code to take advantage of these symmetries is an important component of creating a useful quantum chemistry computational program. A Hartre-Fock Calculation for Water To find the HF MOs we need to solve the HF secular determinant 29-3 (29-1) and find its various roots. We know that we can compute overlap integrals and that Fock matrix elements are defined by F µ! = µ – 1 2 " 2 ! – Z k k nuclei # µ 1 r k ! + P $% $% # µ! $% ( ) – 1 2 µ$ !% ( ) & ’ ( ) (29-2) Finaly, the density matrix elements appearing in eq. 29-2 are defined as P !" = 2 a !i i occupied # a "i (29-3) where the a values are the coeficients of the basis functions in the occupied molecular orbitals. At the first step of an HF calculation, we simply gues what these are, and after that we iterate through solution of the secular determinant to derive new coeficients and we continue until we reach self-consistency. Let's now do this for water. We'll use a structure close to the experimental structure: O−H bond lengths of 0.95 Å and a valence bond angle at oxygen of 104.5 deg. In cartesian coordinates (Å), that puts oxygen at position (0.000, 0.000, 0.116), one hydrogen at (0.000, 0.751, −0.465) and the other hydrogen at (0.000, −0.751, −0.465). Thus, the origin has been taken as the molecular center of mas, the atoms lie in the yz plane (al x coordinates are zero), and the z axis is the symmetry axis of rotation present in the water molecule that makes it belong to the symmetry point group C 2v (don't worry if you don't know anything about symmetry point groups—it won't be important, although it does make the calculation more eficient). The basis functions we wil use are shown in the figure on the next page. Function 1 is the O 1s, function 2 the O 2s, etc. So, with the geometry in hand, we can compute the necesary parts of the secular determinant. We begin with the overlap matrix elements. They are 29-4 H a O H b x y z Basis Functions 1s 2s 2p x 2p y 2p z 1s(H a ) 1 2 3 4 5 6 7 1s(H b ) S = 1.000 0.237 1.000 0.000 0.000 1.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 1.000 0.055 0.479 0.000 0.313 !0.242 1.000 0.055 0.479 0.000 !0.313 !0.242 0.256 1.000 " # $ $ $ $ $ $ % & ’ ’ ’ ’ ’ ’ (29-4) There are many noteworthy features in S. First, it is shown in a lower packed triangular form because every element m,n is the same as the element n,m by symmetry, 29-5 and every diagonal element is 1 because the basis functions are normalized. Note that, again by symmetry, every p orbital on oxygen is orthogonal (overlap = zero) with every s orbital and with each other, but the two s orbitals do overlap (because they are not hydrogenic orbitals—which would indeed be orthogonal as diferent eigenfunctions of a Hermitian operatorthey are just diferent gaussians; S 12 = 0.237). Note that the oxygen 1s orbital overlaps about an order of magnitude les with any hydrogen 1s orbital than does the oxygen 2s orbital, reflecting how much more rapidly the first quantum-level orbital decays compared to the second. Note that by symmetry the oxygen p x cannot overlap with the hydrogen 1s functions (positive overlap below the plane exactly cancels negative overlap above the plane) and that the oxygen p y overlaps with the two hydrogen 1s orbitals equaly in magnitude but with diferent sign because the p orbital has diferent phase at its diferent ends. Finaly, the overlap of the p z is identical with each H 1s because it is not changing which lobe it uses to interact. Now, the kinetic energy matrix is (in a.u.) T = 29.003 !0.168 0.808 0.000 0.000 2.529 0.000 0.000 0.000 2.529 0.000 0.000 0.000 0.000 2.529 !0.002 0.132 0.000 0.229 !0.177 0.760 !0.002 0.132 0.000 !0.229 !0.177 0.009 0.760 " # $ $ $ $ $ $ % & ’ ’ ’ ’ ’ ’ (29-5) Consistent with our analysis from last lecture's homework, every diagonal term is much larger than any off-diagonal term. Note that of-diagonal terms can be negative. That is because there is no real physical meaning to a kinetic energy expectation value involving two diferent orbitals. It is just an integral that appears in the complete secular determinant. Symmetry again keeps p orbitals from mixing with s orbitals or with each other. The nuclear atraction matrix is V = !61.733 !7.447 !10.151 0.000 0.000 !9.926 0.000 0.000 0.000 !10.152 0.019 0.226 0.000 0.000 !10.088 !1.778 !3.920 0.000 !0.228 0.184 !5.867 !1.778 !3.920 0.000 0.228 0.184 !1.652 !5.867 " # $ $ $ $ $ $ % & ’ ’ ’ ’ ’ ’ (29-6) 29-6 Again, consistent with our prior analysis, diagonal elements are bigger than off-diagonal elements. Again, positive values can arise when two diferent functions are involved even though electrons in a single orbital must always be atracted to nuclei and thus diagonal elements must always be negative. Note that the p orbitals al have diferent nuclear atractions. That is because, although they al have the same atraction to the O nucleus, they have diferent amplitudes at the H nuclei. The p x orbital has the smalest amplitude at the H nuclei (zero, since they are in its nodal plane), so it has the smalest nuclear atraction integral. The p z orbital has somewhat smaler amplitude at the H nuclei than the p y orbital because the bond angle is greater than 90 deg (it is 104.5 deg; if it were 90 deg the O−H bonds would bisect the p y and p z orbitals and their amplitudes at the H nuclei would necesarily be the same). Thus, the nuclear atraction integral for the later orbital is slightly smaler than for the former. The sum of the kinetic and nuclear atraction integrals is usualy caled the one- electron or core part of the Fock matrix and abbreviated h (i.e., h = T + V). One then writes F = h + G where F is the Fock matrix, h is the one-electron matrix, and G is the remaining part of the Fock matrix coming from the two-electron four-index integrals. To compute those two-electron integrals, however, we need the density matrix, which itself comes from the occupied MO coeficients. So, we need an initial gues at those coeficients. We can get such a gues many ways—the gues listed below is from a model like Hückel theory that goes beyond just π systems but that is not iterative. It provides: Occupied MO coefficients at cycle 1. 1 2 3 4 5 1 .994311 -.232461 .000000 -.107246 .000000 2 .025513 .833593 .000000 .556639 .000000 3 .000000 .000000 .000000 .000000 1.000000 4 .000000 .000000 .607184 .000000 .000000 5 -.002910 -.140863 .000000 .766551 .000000 6 -.005147 .155621 .444175 -.285923 .000000 7 -.005147 .155621 -.444175 -.285923 .000000 So, MO 1 is very nearly a pure oxygen 1s orbital, MO 2 is mostly a pure 2s orbital with smal contributions from the H 1s orbitals, MO 3 is a bonding combination of the O 2p y and the H 1s orbitals (note that their coeficients change sign because the p orbital also changes sign—that's what makes the overlap bonding!) MO 4 is another bonding orbital betwen the 2s polarized by the 2p z and the H 1s orbitals. Finaly, MO 5 is the pure O 2p x orbital. This is al the occupied orbital coeficients so we can compute the density matrix using eq. 29-3 and we obtain 29-7 P = 2.108 !0.456 2.010 0.000 0.000 2.000 0.000 0.000 0.000 0.737 !0.104 0.618 0.000 0.000 1.215 !0.022 !0.059 0.000 0.539 !0.482 0.606 !0.022 !0.059 0.000 !0.539 !0.482 !0.183 0.606 " # $ $ $ $ $ $ % & ’ ’ ’ ’ ’ ’ (29-7) With P, we can compute the remaining contribution of G to the Fock matrix. We wil not list al 406 two-electron integrals here.. Instead, we wil simply present the Fock matrix as F = !20.236 !5.163 !2.453 0.000 0.000 !0.395 0.000 0.000 0.000 !0.327 0.029 0.130 0.000 0.000 !0.353 !1.216 !1.037 0.000 !0.398 0.372 !0.588 !1.216 !1.037 0.000 0.398 0.372 !0.403 !0.588 " # $ $ $ $ $ $ % & ’ ’ ’ ’ ’ ’ (29-8) So, we're finaly ready to solve the secular determinant, since we have F and S fully formed. When we do that, and then solve for the MO coeficients for each root E, we get new occupied MOs Occupied MO coefficients at cycle 2. 1 2 3 4 5 1 .994123 -.233041 .000000 -.103527 .000000 2 .026666 .834174 .000000 .541326 .000000 3 .000000 .000000 .000000 .000000 1.000000 4 .000000 .000000 .608401 .000000 .000000 5 -.004400 -.131110 .000000 .774771 .000000 6 -.006039 .156851 .443140 -.279656 .000000 7 -.006039 .156851 -.443140 -.279656 .000000 If you compare these coeficients to those above, you can se that they have not changed much in a qualitative way, but there are quantitative changes in many places. So, we iterate again, and again, and again, until we are satisfied that further iterations wil not change either our (i) energy, (i) density matrix, or (ii) MO coeficients (it’s up to the quantum chemist to decide what is considered satisfactory). In our water calculation, if we monitor the energy at each step we find: 29-8 E(RHF) = –74.893 002 803 A.U. after 1 cycles E(RHF) = –74.961 289 145 A.U. after 2 cycles E(RHF) = –74.961 707 247 A.U. after 3 cycles E(RHF) = –74.961 751 946 A.U. after 4 cycles E(RHF) = –74.961 753 962 A.U. after 5 cycles E(RHF) = –74.961 754 063 A.U. after 6 cycles E(RHF) = –74.961 754 063 A.U. after 7 cycles So, our original gues is realy not too bad—off by a bit les than 0.1 a.u. or roughly 60 kcal mol −1 . Our gues energy is too high, as the variational principle guarantes that it must be. Our first iteration through the secular determinant picks up nearly 0.07 a.u., our next iteration an additional 0.000 42 or so, and by the end we are converged to within 1 nanohartre (0.000 000 6 kcal mol −1 ). We can break the energy down into its four components: Electron kinetic energy: 74.606 167 403 Electron-nuclear atraction energy: –197.098 028 755 Electron-electron repulsion energy: 38.265 406 023 Nuclear-nuclear repulsion energy: 9.264 701 265 Note that the first term is and the sum of the next thre terms is . So, we can ask whether the wave function satisfies the virial theorem –2 = . We se that –2 = –149.212 334 806 and = –197.098 028 755 + 38.265 406 023 + 9.264 701 265 = –149.567 921 466. So, the virial theorem is not satisfied exactly (as is typical for a non-exact wave function) but it is satisfied to within about 0.5%. Having converged the SCF equations, we have final MOs from which to form our Slater-determinantal wave function. The final MO coeficients for all MOs (not just the occupied ones) are 29-9 Final MOs: 1 2 3 4 5 EIGENVALUES -- -20.24094 -1.27218 -.62173 -.45392 -.39176 1 .941 -.23251 .000 -.10356 .000 2 .02672 .83085 .000 .53920 .000 3 .000 .000 .000 .000 1.000 4 .000 .000 .6067 .000 .000 5 -.042 -.13216 .000 .7828 .000 6 -.0605 .15919 .4453 -.27494 .000 7 -.0605 .15919 -.4453 -.27494 .000 6 7 EIGENVALUES -- .61293 .75095 1 -.1340 .000 2 .89746 .000 3 .000 .000 4 .000 .9474 5 -.7428 .000 6 -.80246 -.84542 7 -.80246 .84542 where "eigenvalue" refers to the energy of one electron in the MO. Remember that the sum of al of the occupied MO energies should be an underestimation of the total electronic energy because electron-electron repulsion wil have been double counted. So, if we sum the occupied orbital energies (times two, since there are two electrons in each orbital), we get 2(–20.24094 – 1.27218 – 0.62173 – 0.45392 – 0.39176) = –45.961 060. If we now subtract the electron-electron repulsion energy 38.265 406 we get –84.226 466. If we add the nuclear repulsion energy 9.264 701 to this we get a total energy –74.961 765. The diference betwen this and the converged result above (–74.961 754) can be atributed to rounding in the MO energies, which are truncated after 5 places. Notice that the 5 occupied MOs al have negative energies. So, their electrons are bound within the molecule. The unoccupied MOs (caled "virtual" MOs) al have positive energies, meaning that the molecule wil not spontaneously acept an electron from another source. This is certainly in keeping with our expectations. Water neither spontaneously gives up nor grabs electrons. Next time, we wil look more carefully at the orbitals themselves and at the information contained within the wave function. Homework To be solved in clas: Use eq. 29-3 and the initial gues MO coeficients to verify that P 21 in eq. 29-7 is –0.456. 29-10 To be turned in for possible grading December 4: What is the value of P 52 for water after the first SCF step (i.e., at cycle 2)? By what percentage did it change compared to the initial gues? By what percentage did the total energy change for this same step? Laura Gagliardi 3502_29_nove20_LG