J . Phys. Chem. 1988, 92, 733-740
733
Extended Huckel Calculations on Copper Clusters A. T. Amos,* P. A. Brook, and S. A. Moirt Department of Mathematics and Department of Metallurgy and Materials Science, University of Nottingham, University Park, Nottingham NG7 2RD, England (Received: February 4, 1987; In Final Form: August 17, 1987)
A modified version of extended Hiickel theory is used to calculate some properties of copper clusters. The theory seems useful for discussing excitation spectra, spin densities, and ionization potentials of small clusters, but there are difficulties in predicting conformations. The large clusters are intended to model the bulk metal, and the calculated values of their properties are broadly in agreement with the bulk values, although an exception is the cohesive energy where the computed value is much too small.
1. Introduction The purpose of this paper is to report the results of extended Huckel theory (EHT) calculations on copper clusters. Here the term “clusters” is used in two senses-for both calculations on Cuz, Cu,, Cu4, Cu5, tug, i.e., real clusters consisting of a few metal atoms, as well as larger clusters, which are intended to represent the whole of the crystalline metal. The former, especially the dimer, trimer, and pentamer, are of interest because experimental information is becoming available with which to compare and, hence, test theoretical calculations. Large clusters are of importance in the cluster model of surface states and chemisorption. For the cluster model of surfaces to be successful, a large cluster of copper atoms has to possess some of the properties of the bulk metal. This does not appear to be the case for E H T calculations on copper clusters’-3 where the most serious deficiency is that the density of states (DOS) c u y e has far too small a width, even with very large clusters. This is very much the case for those calculations employing a single d orbital but the use of a double { d orbital improves the situation a little. Reviewing this problem, as long ago as 1976, Messmer et aL4 suggested that, to improve this situation, E H T should be reparametrized. Because of the paucity of experimental information, they proposed that the parameters should be chosen so as to reproduce the results of x, calculations, which do appear to be more acceptable. There are no reports in the literature to suggest that this proposal has been taken up. However, if one attempts to do so, it turns out that X, results, of themselves, do not suffice to fix the EHT parameters. Experimental results from Cu2 (which have become available only recently), together with the X, calculations, do enable a better set of EHT parameters to be determined but, even then, the DOS curves for large clusters are still too narrow. It seems, therefore, that, for a real improvement, EHT needs to be made more flexible. Ionization potentials are the counterpart for molecules of DOS for metals and, in view of the failure of EHT to give satisfactory DOS,it is interesting to note that, on the whole, EHT does not produce IP values of sufficient reliability to be useful, for example, in the interpretation of photoelectron spectra. This problem has been addressed by Jerrard et al.s who conclude that the major error is the restricted form of the Wolfsberg-Helmholz approximation for the off-diagonal matrix elements of the E H T Hamiltonian. The use of a more flexible form, together with the empirical determination of the parameters for the matrix elements of the Hamiltonian, improves markedly the IP‘s computed by EHT. In view of this, the same procedure has been adopted for the calculations reported here and, as will be seen, the results for clusters are also much improved. 2. Extended Hbckel Theory and Parameter Determination EHT is a one-electron approximation which considers only the valence electrons of the system. These electrons are assumed to occupy molecular orbitals (MO’s) which are eigenfunctions of an ‘Present address: Corrosion Protection Centre, UMIST, P.O. Box 8 8 , Manchester M60 lQD, England.
effective Hamiltonian defined in terms of its matrix elements with respect to the valence atomic orbitals (AO’s) (w,) by H,,= -I,
H,,= -1.125(Zr,
+ Z’,)SrSexp(-O.l3R)
(1) where R is the distance between the atoms associated with w, and os,respectively, and S,,is the overlap integral between these orbitals. Here we adopt the modification of EHT introduced by Jerrard et al.s whereby different parameters Z’Ji are used for the off-diagonal matrix elements and both sets of parameters are treated as adjustable. The original and commonly used version of E H T has Irr = I , and each I, is set equal to the valence-state IP associated with w,. In addition, we follow Anderson6 by incorporating an exponential term in H,,since this appears to improve calculated bond distances and we also employ the method of Ammeter et al.’ to prevent counterintuitive orbital mixing. In conventional EHT, the total energy E is taken to be the sum of orbital energies, E,, i.e.
E = &,E, (2) where n, is 2, 1, 0, respectively, for doubly occupied, singly occupied, and unoccupied orbitals. Unfortunately, this energy expression does not always have the correct dependence on interatomic distances and angles; for example, in some diatomics, the energy is purely attractive. Anderson6,* has considered this problem and proposes that (2) should be replaced by E = C n i E ,- ZZbJp,(r)lRb
- rl-’ dr
(3)
where the sum is over all pairs of atoms a and b (taking a to be the more electronegative) and Zb is the nuclear charge of atom b, Rb its position, and pa the charge density of the free atom a as determined from its valence orbitals with an equivalent nuclear contribution. The additional term in (3) represents the repulsive interaction which arises from the superposition of rigid atoms and its inclusion generally improves calculated potential energy curves. We have used (3) for the present calculations. As was pointed out in the Introduction, Messmer et a].: when discussing deficiencies in EHT, suggested that more satisfactory results might be obtained by parametrizing E H T calculations against X, calculations, since the latter appeared to provide the most successful treatment of cluster models. We have adopted this excellent suggestion but, in practice, we have found it necessary to use in addition experimental information on the copper dimer. The parameters thus obtained are listed as parameters A in Table ~~~
(1) Baetzold, R. C. J. Phys. Chem. 1978, 82, 738. (2) Baetzold, R. C.; Mack, R. E. J . Chem. Phys. 1975, 62, 1513. (3) Anderson, A. B. J. Chem. Phys. 1978, 68, 1744; Transition Metals; Conf. Ser. No. 39; Institute of Physics: Bristol, U.K., 1978; p 379. (4) Messmer, R. P.; Knudson, S. K.; Johnson, K. H.; Diamond, J. B.; Yang, C. Y.Phys. Rev. B 1976,813, 1396. ( 5 ) Jerrard, R. J.; Amos, A. T.; Brook, P. A. Theor. Chim. Acta 1982,61, 163. ( 6 ) Anderson, A. B. J . Chem. Phys. 1975, 62, 1187. (7) Ammeter, J. H.; Burgi, H. B.; Thibeault, J. C.; Hoffmann, R. J. A m . Chem. SOC.1978, 100, 3686. (8) Anderson, A. B. J. Chem. Phys. 1974, 60, 2477.
0022-365418812092-0733$01.50/0 0 1988 American Chemical Society
734 The Journal of Physical Chemistry, Vol. 92, No. 3, 1988 TABLE I: Extended Hiickel Parameters for CoppeP parameters A I ; , = -6.4 I ; , = -3.6
14s = -7.0 14, = -3.1
Amos et al.
parameters B
I,, = -6.9 14, = -4.9
parameters A and B
P4,= -6.0 '4,
I,, = -8.5
= -5.0
{dS
I;, = -14.5
= 1.75
{4:4p
= 1.80
+
OThe d orbital is 0 . 5 9 3 3 2 ~ ~0.57442wb, where w, and wb are 3d atomic orbitals with exponents 5.95 and 2.3, respectively.
TABLE 11: Experimental and EHT Calculated Values for Various Properties of Cu2 calcd values property
exptl value
bond length/A binding energy/eV ionizn potential/eV stretch freq/cm-'
2.22" 2.04 i 0.13; 2.0Y 6.7c 7.37: 7.89e 266"
HF' 2.41 0.68 6.37 1979
CE' 2.40 2.07
200
x,g
2.22 2.10 7.35h 286
EHT(A)' 2.26 2.04 7.83 317
EHT(BY 2.20 2.29 7.64 348
OReference 13. bReference 14. 'Reference 15. dReference 16. eReference 17. /Reference 18. #Reference 19. h T h e calculated value was 4.35 but this was for a ground-state calculation. A transition-state calculation would increase it, probably by about 3 eV as shown in the table. 'Parameters A. 'Parameters B
I and we now discuss some of the more important considerations which led to this particular choice. When these parameters are obtained, an important factor has been the equilibrium bond length & of the dimer. The parameter which has the major effect on the EHT predicted value for Ro is the exponent c4, of the 4s valence atomic orbital but rather large values of this are needed to reproduce the experimental result. Good quality Hartree-Fock and even CI calculations, both of which usually give accurate bond lengths, overestimate & for Cu2 by about 0.04 A. This has led to detailed investigations by several authors (the most recent, which has references to earlier work, is ref 9) and the conclusion reached is that relativistic effects are responsible for a shortening of the Cu2bond by about this amount. The implication is that, where relativistic effects are not included, as, of course, will be the case in EHT, a good theoretical prediction would be Ro = 2.26 A, rather than the experimental Ro = 2.22 A. Bearing this in mind, therefore, we set {4s at 1.75, which does cause EHT to predict the slightly longer length. An alternative point of view is that a semiempirical theory should be allowed the complete flexibility to predict properties such as Ro by means of a suitable choice of parameters which would take into account, through that choice, relativistic effects on the core potential. Normally we would entirely accept this point of view. That we dissent in this particular case is because we feel the dimer equilibrium distance is somewhat anomalous and fitting to it would cause bond lengths in the copper clusters to be underestimated. Experimental bond lengths for small copper clusters would help resolve this point. The measured IP of the copper dimer can be of help in finding suitable diagonal and off-diagonal 4s matrix elements since the IP will be -Eiwhere E, is the energy of the highest occupied MO. This orbital is almost certainly of crg symmetry, built up mainly of 4s atomic orbitals, so that its energy will depend on Z4, and I;,. Although the IP alone will not fix both of these parameters, there is an excited-state u,, energy level, essentially the antibonding associate of the ug orbital, which is also built up mainly of 4s atomic orbitals. If its energy were known, it would be possible to determine both Z,, and I;,. An X, calculation on Cu2 by Ozin et al.,IO and also Hartree-Fock calculations," suggest that the energy difference between the ugand the u,,orbital might be about 2.4 eV. Taking this value, the Z4s and P4, parameters (after rounding to two significant figures) will be those given in Table I. The energies of the unoccupied orbitals of mainly 4p character are not known with any precision. It has been suggested'* that
the lowest of these has an energy below the uu 4s-type orbital, but this does not seem to be generally accepted. Until this point has been definitely settled, perhaps by an unambiguous interpretation of the excitation spectrum, it seems best to proceed by assuming that the 4p orbitals have energies above the uu orbital, implying that Z4p = = -3 to -4 eV, as listed in Table I for parameters A. For the 3d valence atomic orbital, we have used a double { orbital, which is generally regarded as more satisfactory than a single Slater orbital for representing d atomic orbitals. The values for the 3d parameters were determined mainly by a least-squares fit to the X, calculations4 on Cu13.Significantly, we found that in order to reproduce the d-band width it was essential to take lZ'3dl to be much greater than lZ3dl (see Table I), thus increasing the absolute values of the off-diagonal elements H,, involving d orbitals. A referee has pointed out that E H T calculations with two d orbitals give broader bands than those with one. This is an interesting point since it suggests that an alternative way to increase the values of the off-diagonal elements would be to increase the overlap integrals by including a very diffuse second d orbital. Table I1 shows the agreement between the EHT values, denoted by EHT(A), obtained for the copper dimer by using parameters A of Table I and experiment. For comparison, some other calculated values are included also. Of course, the EHT(A) results do agree very well with the experimental ones, because the parameters were chosen with this in view. The slight discrepancies which remain are due to rounding the parameters at two significant figures. Note, however, that no attempt was made to fit the fundamental vibration frequency, so it is pleasing that the EHT(A) value is the correct order of magnitude. The discrepancy of about 20% is typical of that found with calculations based on singledeterminantal wave functions. In section 4 we show that parameters A lead to what are probably incorrect predictions of the conformations of small copper clusters. A referee has suggested that we should exapine whether this can be improved by bringing the 4p level closer to the 4s level so that the 4p orbitals contribute more to the occupied u-bonding orbitals. Unfortunately, this cannot be done in isolation or else binding energies become too large. To counteract that the 4s parameters, particularly Z;,, can be changed, although this tends to lead to poor IP values. A compromise is represented by the parameter set f3 listed in Table I. Note that the orbital exponents and d-orbital parameters have been left unchanged for simplicity. The dimer properties calculated with this set are given in Table
(9) Scharf, P.; Brode, S.; Ahlrichs, R. Chem. Phys. Lett. 1985, 113,447. (10) Ozin, G. A.; Huber, H.; McIntosh, D.; Mitchell, S.; Norman Jr., J. G.; Noodleman, L. J . Am. Cfiem. SOC.1979, 101, 3504. (1 1) Miyoshi, E.; Tatewaki, H.; Nakamura, T. J . Chem. Phys. 1983, 78, 815. (12) Kolb, D. M.; Rotermund, H. H.; Schrittenlacher, W.; Schroeder, W. J . Cfiem. Pfiys. 1984, 80, 695. (13) Aslund, N.; Barrow, R. F.; Richards, W. G.; Travis, D. N. Ark. Fys.
(14) Gingerich, K. A. Cum. Top. Mat. Sci. 1980, 6, 345. (15) Kant, A. J. Chem. Pfiys. 1964, 41, 1872. (16) Cabaub, B. Ph.D. Thesis, University of Claude Bernard, Lyon, 1972, as quoted in Joyes, P.; Leleyer, M. J . Phys. E 1973, 6, 150. (17) Powers, D. E.; Hansen, S. G.; Geusic, M. E.; Michalopoulos, D. L.; Smalley, R. E. J. Cfiem. Phys. 1983, 78, 2866. (18) Shim, I.; Gingerich, K. A. J . Chem. Pfiys. 1983, 79, 2903. (19) Delley, B.; Ellis, D. E.; Freeman, A. J.; Baerends, E.-J.; Post, D. Pfiys. Reo., B 1983, 827, 2132.
1965, 30, 171.
Extended Huckel Calculations on Copper Clusters 11, denoted by EHT(B). They are in reasonable agreement with experiment and this could probably be improved by a further refinement of the parameter values. However, the overestimate of the binding energy and underestimate of the IP are so interrelated that it would be difficult to improve one without making the other much worse. These difficulties would be accentuated if the 4p level were reduced further. As it is, it has been lowered just enough to give a good prediction of the trimer geometry. Usually it is assumed that a single EHT calculation will suffice to determine a multitude of a molecule's properties, including ground-state properties such as energy (and, hence, conformations) and charge distributions; excitation energies involving transitions to excited states; IP's and DOS arising from transitions into final states with one electron fewer than the ground state. A theory based on a initial-state wave function cannot be expected to do all of this equally well. For example, in spite of Koopman's theorem, orbital energies from Hartree-Fock calculations on copper clusters do not correlate at all with IP's or with bulk DOS (the Hartree-Fock d band is well-separated from the s band). Thus, it would seem that, if E H T is used to predict IP's and if its parameters are based on experimental IP's, it should not be thought of as an initial-state theory and certainly not, in any sense, an approximation to Hartree-Fock theory. For excitation-energy calculations, Mehrotra and HoffmannZ0have suggested that EHT should be interpreted as a transition-state theory for transitions between states with an equal number of electrons. Following this idea, Jerrard et aLs proposed that for E H T IP calculations the transition state should be that for transitions between states whose numbers of electrons differ by one. If this is correct then parameters A are clearly applicable to such a transition-state Hamiltonian but it could well be the case that the orbital energies from this Hamiltonian are not the proper ones to use in the energy formula (3) and that parameters similar to set B defines a more suitable ground-state Hamiltonian with more appropriate orbital energies.
3. The Excitation Spectrum of Cu2 In the simplest form of EHT, excitation energies are given as the differences between orbital energies. As just stated, it has been arguedZothat, for this to be correct, the Hamiltonian would have to correspond to a transition-state Hamiltonian for transitions between states with the same number of electrons which would differ somewhat from the Hamiltonian whose energies are related to I P S . This is precisely what is found in the X, method, where a different Hamiltonian is used to discuss excitation energies compared with that used for ionization energies. Assuming the same should apply in EHT, the best approach would be to find parameters to define the appropriate Hamiltonian for the excitation spectra of copper clusters but not sufficient of the experimental information which would be needed to obtain these new parameters is available. These circumstances mean that the old parameters must be used, even though they are not entirely suitable, with the consequence that qualitatively correct results are all that can be expected. Nevertheless, it seems worthwhile proceeding along these lines, since the interpretation of the spectra of small metal clusters is difficult, and even results that are only qualitatively correct can help. The EHT(A) orbital energies and orbital symmetries for the dimer are listed in Table I11 and, for comparison, the orbital energies obtained by Ozin et a1.I0 using the X, method are also shown. The X, results are specifically for excitation-energy calculations and not for IP's and, therefore, have to be decreased by about 3.5 eV for comparison with the EHT values. When this is done, the two sets of results are in very reasonable agreement. Moreover, the significant factor is not the absolute vaues, but the relative values, since the excitation energies will be given by the differences between the energies of occupied and unoccupied orbitals. Therefore, the excitation spectrum predicted by Ozin et al. is compared in Table IV with the two sets of EHT results, obtained by using the two parameter sets. (20) Mehrotra, P. K.; Hoffmann, R. Theor. Chim. Acfa 1978, 48, 301.
The Journal of Physical Chemistry, Vol. 92, No. 3, 1988 735 TABLE III: EHT and X, Orbital Energies and Symmetries for Cu2 EHT results (Darameters A) s, P, d component % ~
energy/eV sym unocc orbitals
occ orbitals
+12.80 -2.10" -3.37 -3.52" -5.57 -7.83 -8.11" -8.12 -8.42" -8.58" -8.97" -9.35
3uu 28 2a, 0 3ug 1 2 ~ , 0 2u, 68 2ug 87 0 lk, la, 4 0 16, 16, 0 la, 0 la, 12
X, resultsb
p
d
energy/eV
sym
68 98 98 99 31 0 2 1 0
4 2 1 1 1 13 98 94 100 100 99 86
0.20 -0.84" -1.99 -4.38 -4.42" -4.74 -4.77" -5.00" -5.42" -6.48
3a, 2a, 20, la, lr, 2a, 16, 16, In, la,
s
0 1 2
Indicates doubly degenerate orbital. Reference 10. TABLE I V Theoretical Spectrum of Cu2
theor excitn energy/eV EHT(A)" 2.26 2.54 3.78 4.31 4.59 5.06 4.60 2.55 3.01 2.85 3.40
EHT( B ) ~ 1.52 2.02 3.18 2.24 2.74 3.18 2.85 2.01 2.46 2.29 2.90
XeC 2.75 2.43 4.49 3.90 3.58 4.16 4.58 2.39 3.01 2.78 3.43
excited state 2ug
--
2a,(Z,) 2o"(n,) la, 2U" 2ag 2a, 2a, l a , -.+ 16,- 27ru 3 5 1U" la, 2U"d 2afi 16, 16, 2U"d Is, 2ufi lag
+
----
"Parameters A. Parameters B. Reference 10. Dipole-forbidden transitions.
On the whole, the agreement between EHT(A) and the X, results is quite good, but there are systematic differences. The energy of the 2?r, orbital in the X, calculation is much closer to the energy of the occupied orbitals than is the case in EHT(A) and this accounts for the discrepancies in lines 4-6 of Table IV. These two calculations predict a different ordering for the lowest transitions as the 2?r, orbital is predicted by EHT(A) to be the highest occupied orbital but it lies among the d orbitals in the X, calculation. Experimentally, two bands have been observed at 2.53 and 2.69 eV which are probably transitions to the 8,and nuexcited states but the ordering has not been determined experimentally as yet. However, Aslund et al.13 claim that an analysis of the vibrational-rotational spectrum indicates that the upper level is 2, but, on the other hand, the force constants indicate that the energy curve for the excited state is very flat, which would suggest it is the II, state. A CI calculation by Witko and Beckman# predicts the lowest excited state to be 8,at an energy of 2.91 eV above the ground state with the Z, level at an energy of 0.38 eV above the u state. Unfortunately, the accuracy of the calculation is not sufficient to be absolutely certain that this ordering is correct. The EHT(B) values for the excitation energies are much smaller than those found by using set A, largely because the EHT(B) orbital energies for the unoccupied orbitals (which are mainly of p character) are lower in energy. The lack of agreement with experiment shows that parameter set B is not suitable for calculating excitation energies. 4. Geometry and Properties of Small Copper Clusters A topic of recent theoretical interest has been the geometry of small copper clusters containing up to six Cu atoms. There have been a number of EHT calculations which have found that, in most cases, the linear structure is the most favored on energy (21) Witko, M.; Beckmann, H. 0. Mol. Phys. 1982, 47, 545.
136
Amos et al.
The Journal of Physical Chemistry, Vol. 92, No. 3, 1988 IPleV
TABLE V: Atomization Energies per Atom of Cu, ( n = 2-6) Structures As Predicted by Various Theoretical Methods (in eV) n structure' EHT(A)b EHF(B)' HFd DIM'
D-=h C2h
1.02 (1.02) 0.78 0.91 0.94 (0.98) 0.66 0.94 1.04 1.07 1.11 0.88 0.89 0.88 1.01 1.08 0.92 1.07 1.15
1.14 (1.02) 1.21 1.17 1.22 (0.98) 1.59 1.40 1.28 1.50 1.50 1.23 1.61 1.52 1.30 1.58 1.62 1.34 1.60
0.42 (1.02) 0.39 0.43 0.39 (0.98) 0.44 0.52 0.55 0.50 0.50 0.60 0.61 0.54
75
0.99 (1.02) 1.07 1.01 1.09 (0.98) 1.60 1.56 1.20 1.58 1.22 1.22 1.76 1.78 1.20
55 40
PY.
75
60
80
1cc)
120
140
160
180
%-age spin density
toot
A
x at central atom o at end atom
2.17 1.27 1.30
a C,, (n odd) and c*,+ (n even) are planar zig-zag structures. D3* is an equilateral triangle for n = 3 and a trigonal bipyramid for n = 5. D-,, is linear. DZhis a rhombus. Ddhis a square when n = 4 and a body-centered square when n = 5. rh is a tetrahedron. C , is a square pyramid. Oh is an octahedron. *Parameters A. 'Parameters B. dReference 22. eReference 23. fReference 14. gQuoted in ref 24.
ground^.^ Demuynck et aLz2 have used Hartree-Fock wave functions and predict a change from a linear to a three-dimensional structure at Cu5. Richtsmeler et al.23using the method of diatomics in molecules (DIM), based on Heitler-London theory, find that the three-dimensional structures have the lowest energies. Both sets of EHT parameters from Table I have been used to compute atomization energies, IPS,and equilibrium configurations for various structures and the results are given in Tables V and VI. In Table V, the atomization energies are compared with experimental and other theoretical values, where available, and it is clear that EHT(A) results differ considerably from those of the DIM method. The EHT(A) calculations favor two-dimensional zig-zag structures as compared to the compact structures which are predicted to be the most stable in DIM theory. Moreover, the EHT(A) atomization energies are generally lower than the DIM energies and the Hartree-Fock energies are smaller still, as is to be expected. On the other hand, when parameters B are used, there is greatly improved agreement with the DIM values. Moreover, the three-dimensional structures are now favored. Other experimental information is available for the trimer and pentamer from Howard and collaborators who have measured the ESR spectrum of Cu, (ref 25) and Cu, (ref 26) trapped in matrices. Howard, Preston, et al. interpret their ESR results for Cu3 as indicating 4s unpaired spin populations of 29% for each of the terminal copper atoms and 3% for the central copper atom. Because of this, they claim the molecule must be slightly bent. In D3hgeometry (an equilateral triangle) the highest occupied energy level, which is the energy of the unpaired electron, is twofold degenerate. The degenerate orbitals are mostly built up of 4s atomic orbitals with only a small admixture of d atomic orbitals and, when parameters A are used, they have the following form in terms of the 4s orbitals wl, w 2 , and w 3 : 4, = 0 . 6 2 ~- ~0 . 6 2 ~ ~ 42 = 0 . 3 2 ~ 1- 0 . 7 2 ~ 2+ C . 3 6 ~ 3
I
(4)
As the bond angle is changed from 0 = 60°, so that the molecule takes up C, (isosceles triangle) geometry, this degeneracy is split. For 0 > 60°, the 4, orbital is of lower energy and thus is the orbital (22) Bachmann, C.; Demuynck, J.; Veillard, A. Gazz. Chim. Italiana 1978, 108, 389; Faraday S y m p . Chem. SOC.1980, 14, 170. (23) Richtsmeler, S. C.; Hendewerk, M. L.; Dixon, D. A,; Gole, J. L. J , Phys. Chem. 1982, 86, 3932.
Figure 1. Variation of ionization potentials and spin densities with angle in the copper trimer.
occupied by the unpaired electron, while for 0 < 60' it is 42which is lower and therefore becomes occupied. The forms of these orbitals remain much the same when the angle is changed and, therefore, the spin density will be very different, depending on whether the unpaired electron occupies orbital q51 (60' < 0 < 180') or orbital 42 (Oo < 0 < 60'). The spin densities of the central and end atoms are plotted as a function of B in Figure 1. Obviously, it is only for 0 > 60' that there is agreement with the ESR spectrum. However, this agreement holds for any value of 0 greater than 60' and not only, as assumed by Howard, Preston, et al., for values near 180'. The DIM method predicts two equilibrium structures for Cu,, which are Jahn-Teller distortions of the equilateral-triangle configuration and have angles 0 = 65.6' and 0 = 55.4'. The Hartree-Fock calculations of Bachmann et aLz2as well as Anderson's EHT result2' predict linear geometry. More recent Hartree-Fock results, however, agree with DIM in finding the lowest energy structure is an isosceles triangle with 0 = 77.6' l i and 0 = 53.9°.24 As Table VI shows, parameter set B produces a result in broad agreement with this, predicting an angle of 0 = 62O but, of course, the parameters were chosen specifically to do this. On the other hand, parameters A predict a much larger angle of 0 = 124', which is almost certainly incorrect. Since the major difference. between these parameters is the 4p level, it seems likely that polarization functions play a crucial role in determining the equilibrium geometry. In their work on supersonic copper clusters in the 1-29-atom range, Powers et al.17were able to obtain information on the I P S by using a number of fixed frequency lasers to examine the threshold behavior of the photoionization cross section. They conclude that all clusters have an IP of less than 7.9 eV and greater than 4.98 eV. In a specific study of Cu3 by the same its IP is found to be 5.48 f 0.5 eV. Figure 1 shows the EHT(A) value of the IP as a function of the bond angle 8. Although the experimental bounds are large, it is clear that it is only for angles near 60° that the EHT(A) estimates agree with experiment. I t is clear also from Table VI that the IP predicted with parameters B is too large. To summarize the situation with Cu3, therefore, the present EHT results agree with DIM and most recent Hartree-Fock (24) del Conde, G.; Bagus, P. S . ; Novaro, 0. Phys. Rev. B 1982, B25, 7843. (25) Howard, J. A,; Preston, K. F.; Sutcliffe, R.; Mile, B. J . Phys. Chem. 1983, 87, 536. (26) Howard, J. A,; Sutcliffe, R.; Tse, J. S.; Mile, B. Chem. Phys. Lett. 1983, 94, 561. (27) Morse, M. D.; Hopkins, J. B.; Langridge-Smith, P. R. R.; Smalley, R. E . J . Chem. Phys. 1983, 79, 5316.
The Journal of Physical Chemistry, Vol. 92, No. 3, 1988 137
Extended Hiickel Calculations on Copper Clusters
TABLE VI: EHT Geometries and Ionization Potentials for Various Structures of Cu. ( n = 3-6) EHT parameters A
structure'
n 3
bond length/A 2.41 2.33 2.32 2.48 2.37 2.34 2.41 2.3 1 2.36 2.38 2.49 2.35 2.32 2.50 2.35 2.32
D3h D-h
CZ" Td D4h D-h D2h CZh
5
D4h C4" D3h
D- h C2" 6
bond angle/deg
O h
D- h
c?.h
EHF parameters B bond lenath/A bond angle/deg
IP/eV 5.77 7.12 6.87 5.83 6.48 7.56 6.80 7.37 7.03 6.33 5.8 1 7.13 6.79 6.41 7.45 7.16
124
120 121
122 120b
2.30 2.24 2.28 2.36 2.28 2.25 2.32 2.32 2.27 2.24 2.38 2.25 2.33 2.39 2.25 2.34
62
1 206 60b
60b 60b
IP/eV 6.33 7.28 6.35 6.34 6.79 7.50 7.03 6.39 7.12 6.68 6.43 7.30 6.73 6.79 7.45 6.76
"See footnote a to Table V. bAssumed value, not optimized. TABLE VII: Variation of Bond Length with Cluster Size Using Parameters A
no. of atoms
Distorted CLV
Distorted D3,,
2 3 4
Distorted 04h
5
.
32
/
26
22
26
0
34
0
.
32
6 14 23
cluster shape linear equilateral triangle tetrahedron bipyramid octahedron copper unit cell see Figure 4
equilib band length/A 2.26 2.40 2.44 2.49 2.50 2.55 2.57
+
Dmh
Figure 2. Spin density distribution in various conformations of the copper pentamer.
calculations in predicting a bent structure. However, the bond and angle depends critically on the values used for 14pand only the parameters B give what is probably the correct angle, Le., near to 60'. Predicted spin densities are insensitive to this angle and to the parameter set, and qualitative agreement with the measured ESR spectrum is found for any angle greater than 60'. The I P is known to be in the range 5.0-6.0 eV and this supports a bent structure with an angle near 60°, but only parameters A give an I P value in the correct range. On the basis of this, it would appear difficult to find a single set of EHT parameters which simultaneously give satisfactory geometries and IPS. The ESR spectrum for Cus has been interpreted as due to a pair of equivalent atoms with about 30% unpaired 4s spin population on each. The remaining three atoms, two of which are equivalent, have very little unpaired spin density. Howard, Sutcliffe, et a1.26 suggest that this implies the structure of Cus is a distorted trigonal bipyramid, Le., the C , structure. To attempt to confirm this, the 4s spin densities for the structures in Table V have been computed by using parameters A and are shown in Figure 2. Since Jahn-Teller distortion is likely to prevent the C,, D3h,and D4hstructures from being totally symmetric, the spin densities at a slightly distorted conformation are the ones given in Figure 2. I t is clear that any of these distorted structures have a spin density distribution in accord with experiment, while the linear and zig-zag structures do not. When IP's are considered, the work of Powers et al." implies that, at least for the small clusters Cu2 to Cus, those with an even number of atoms have an ionization threshold above 6.4 eV and those with an odd number of atoms an ionization threshold below this. The EHT values for IPScorresponding to various structures are listed in Table VI. The IP's for Cu5 found by using parameters B are all above 6.4 eV, in conflict with experiment. This is a further indication that it is set A which is more suitable for predicting IP's. Considering the EHT(A) results only, we note
that for c u 4 and cu6 the majority of values are above 6.4 eV, a notable exception being the tetrahedral structure for Cu4 which implies that this cannot be the most stable configuration. Even for Cu6 the octahedral structure is only just above 6.4 eV and probably this too cannot be the preferred configuration. For Cu5 only the D3h and C4, structures have IP's below 6.4 eV and, therefore, this is some evidence that one of these is the most stable conformation. The tentative conclusion, therefore, is that odd-atom clusters prefer the highly symmetric three-dimensional structures, whereas the even-atom ones take up the more open two-dimensional ones. However, it should be noted that these conclusions are not always in accord with the values for the atomization energies in Table V.
5. Large Copper Clusters The parameters are now used for calculations on large copper clusters, beginning with clusters containing relatively few copper atoms and increasing to 63, the effective limit for the computer available in Nottingham. Except where stated, parameters A have been used. The variation of Cu-Cu bond length with cluster size was considered first and Table VI1 lists the optimized bond length for various clusters, beginning with very small highly symmetric structures and continuing up to a 23-atom cluster. Unfortunately, it is a lengthy process to obtain optimum geometries and it was not practicable for us to optimize the bond length for clusters larger than this. Table VI1 shows an initial large increase in the equilibrium bond length, but between the 6- and 23-atom cluster there is only a small change, suggesting that the bond length has almost converged. The value of 2.57 8, for the largest cluster considered is almost identical with the measured value for the bulk which is a pleasing result. Although the results in Table VI1 are consistent with each other and it is not unreasonable to present them in this light, it is only fair to state that it is easy to find clusters where the equilibrium bond length is very different from those given in the table. For example, an eight-atom cluster in the shape of a cube was found to have a bond length of only 2.41 A. Howcver, it can bc argued that this cluster, and others like it, do not give a good representation of the conformation of the copper atoms in a metal and
738 The Journal of Physical Chemistry, Vol. 92, No. 3, 1988
Amos et al.
TABLE VIII: EHT Calculated Properties (in eV) of Some Copper Clusters no. of atomizn energy energy of d-band Cu atoms per atom HOMO width" EF - Ed,6 19(A)' 23(A)' 41(A)' 63(A)' 23(B)d 41(B)d 63(B)d
bulk values
1.22 1.23 1.26 1.28 1.88 1.98 2.14 3.5'
5.93 5.95 5.63 5.72 6.58 6.45 6.44 4.46
2.36 2.42 2.48 2.66 2.35 2.28 2.48 3.15
1.52 1.55 1.86 1.80
1.07 1.14 1.08 1.6
"The d band is taken to consist of all orbitals with more than 80% d composition. b E Fis the energy of HOMO and Ed" is the upper edge of the d band. 'Calculated by using parameters A. dCalculated by using parameters B. 'The cohesive energy of the metal. 'The experimental work function. that only clusters which do so ought to be considered. Of course, this criterion is, to some extent, subjective but, on the whole, it is easy to identify unsuitable clusters and the problem tends to disappear as the size of the cluster increases. In the remainder of this paper, the results are obtained with the nearest-neighbor bond length fixed at the bulk value of 2.56 A. With this bond length, therefore, and using quite large clusters which, it is hoped, give a good representation of the metal, the calculated properties given in Table VI11 were obtained and they are compared with bulk value^.^^*^^ In making this comparison, it should be recognized that many of the properties will be influenced by final-state effects such as shake-up and plasma excitations. As indicated in the final paragraph of section 2, we believe that EHT, when based on parameters obtained from IP data (e.g., parameters A), has to be thought of as a transition-state theory which, therefore, has final-state effects built into it. If this were not the case, a comparison with bulk values would be meaningless. Considering first the results obtained by using parameters A, we see that a very unsatisfactory feature of this table is the atomization energy per atom which appears to converge to a value of 1.28 eV which is considerably below the cohesive energy of 3.5 eV per atom for the bulk. To some extent this is an unfair comparison since the clusters contain fewer pairwise bonding interactions than does the solid. It has been suggested30 that to allow for this in cluster calculations the cluster atomization energy should be scaled upwards by a factor which is the ratio of the average coordination number of the atoms in the solid to that of the atoms in the cluster. Applying this scaling would improve our estimate of the bulk cohesive energy but it would still be much lower than the true value. The scaling procedure can be regarded as an attempt to extrapolate to infinitely large clusters but it seems unlikely that the size of the cluster is the critical factor. The energy of the highest occupied MO(HOMO), which corresponds to the Fermi level, decreases with increasing cluster size and the width of the d band increases. It is disappointing that neither seems to have converged, even with the very large clusters considered here. Nevertheless, it is encouraging that, even unconverged, the width of the d band is much closer to the experimental value for the bulk and very much larger than is found in more conventional EHT calculations. In Table VI11 it has been assumed that the experimental equivalent of the energy of the HOMO is the work function. If this is correct, it can be seen that the calculated energy is about 1 eV too large, perhaps implying that a further refinement of the E H T parameters is needed. More than any other property, the energy difference between the Fermi level and the onset of the d band seems to depend on the configuration, rather than the size of the cluster. The correct (28) Kittel, C. Introduction to Solid Stare Physics, 4th ed.; Wiley: New York, 1971. (29) Moruzzi, V. L.; Janak, J. F.; Williams, A. R. Calculated Electronic Properties of Metals; Pergamon: New York, 1978. (30) Bicerano, J . ; Keem, J. E.; Schlegel, H. B. Theor. Chim.Acta 1986, 70, 265.
,Ev-10
-9
-8
-7
-6
-5
I
Figure 3. Density of states curve for a 23-atom cluster in the form of a double cube (each cube corresponding to a unit cell with Cu atoms at the corners and centers of each face). The dotted vertical line is at the Fermi
level. value is thought to be about 1.6 eV which is somewhat smaller than that obtained for the larger clusters. It must be understood, however, that in EHT, where there is overlap between the AO's, many of the MO's are a mixture of s, p, and d orbitals and, therefore, it is not easy to decide precisely where the upper edge of the d band occurs. Turning now to the results from parameters B, we find that for most of the properties in Table VI11 (in effect, those related these results are much less satisfactory than the EHT(A) to IPS) ones. A notable exception, however, is the atomization energy which is appreciably larger than the EHT(A) value. Taking into account the fact that it has not completely converged for the 63-atom cluster, it is probable that the EHT(B) estimate of the cohesive energy would be about 1 eV larger than the EHT(A) estimate and much closer to the experimental value although still too small by about 1 eV. The comparison of the values in Table VI11 with bulk values is assuming that the fully converged values will not differ qualitatively from those found for the larger clusters. Unfortunately, there is no comprehensive theoretical treatment which examines the convergence toward bulk values with increasing cluster size for the various computational methods. However, there have been a few calculations from which some tentative conclusions can be drawn. Baetzold' has compared computed EHT DOS curves for Cu clusters containing up to 43 atoms with the bulk DOS computed with the same parameters. H e finds that the width of the bulk d band is about greater than that of the clusters although with his parameters the bulk value is only 2 eV. Being optimistic and applying the same factor to the present results would give a bandwidth in reasonable agreement with experiment. Baetzold notes a distinct variation of the HOMO energy with cluster size but the variation in our calculations appears to depend, also, on cluster shape. Davidson and Fain3' have made calculations which show that changes in the IPSof clusters can still be found on progressing from 79 to 135 atoms. This suggests that extrapolation to work function and Fermi level values from cluster HOMO energy values is difficult and, therefore, the comparisons we have made with experimental bulk values for these quantities must be regarded as somewhat tentative. The discrete energy levels obtained in cluster calculations can be broadened so as to produce DOS curves. To do this we have used a Gaussian function of the form c
where the broadening parameter u was set equal to 0.1 eV. The resulting DOS curves obtained for the 23-, 41-, and 63-atom (31) Davidson, E. R.; Fain, S . C. J . Vac. Sci. Technol. 1976, 13, 209.
Extended Huckel Calculations on Copper Clusters
Figure 4. Density of states curve for a 41-atom cluster in the form four cubes in line, each cube being a unit cell.
,Ev-ll
-10
-9
-8
-7
-6
The Journal of Physical Chemistry, Vol. 92, No. 3, 1988 739
Figure 6. Density of states for the d, s, and p band in the 63-atom cluster. Full line is the d band, the dotted line is the s band. and the broken line is the p band.
-5
I
Figure 5. Density of states curve for a 63-atom cluster in the form of eight cubes, in a 2*2*2 configuration, each cube being a unit cell.
clusters are shown in Figures 3-5. These clusters are made up of two-, four-, and eight-unit cells, respectively. The four unit cells of the 41-atom cluster are in a straight line while the eight-unit cells of the largest cluster are in a symmetrical 2*2*2 array. The DOS curves can be compared with the experimental photoelectron spectrum (see, for example, ref 32). Because the photoelectron spectrum does not have a high resolution, it is only possible to say that the overall shape of the EHT DOS curves is correct. A more telling comparison is with DOS curves from band-theory calculations (see, for example, refs 33, 34) since these display more structure. Although the EHT curves are rather more peaked, we would claim that, on the whole, there is reasonable agreement with the band-theory DOS curves. It can be useful to break down the total DOS into various components such as the d-band, s-band, and p-band density of states. In order to do so the percentage of s, p, and d occupancy of each MO must be determined and these are not easily defined in the case of nonorthogonal AOs. The choice of definition finally adopted is described in the Appendix. The results in Figure 6 show, as expected, a relatively narrow and peaked d band with a much wider and flatter s band overlying it. The p band, of course, lies at much higher energies. Even a cluster of 63 atoms has most of its atoms on the outside of the cluster and so the cluster DOS curves represent mainly the surface DOS. Arlinghaus et al.35have made fully self-consistent
Figure 7. Density of states of surface and interior atoms in the 63-atom cluster. The narrower curve is due to the surface atoms; the wider curve is due to the interior atoms.
surface electronic structure calculations which demonstrate that the surface DOS is narrower than that of bulk and is pushed closer to the Fermi energy. Slab calculations by Citrin et al.36 have reproduced this effect and analyses of valence band XPS spectra appear to confirm the theoretical r e s ~ l t . It~ is~ interesting, ~ ~ ~ therefore, to divide the EHT cluster DOS into two components, that associated with those atoms on the outside of the cluster, which should correspond to surface atoms, and that associated with interior atoms, corresponding to the bulk. This has been done for the 63-atom cluster and the results are displayed in Figure 7 . The surface DOS curve is narrower and shifted toward the Fermi level and, therefore is entirely in accord with the selfconsistent calculations. The interior DOS curve overall seems to be in rather better qualitative agreement with the band-structure DOS curve than is the DOS curve for the full cluster. This demonstrates the importance of separating the total cluster DOS into two components for “bulk”-like and “surface”-like atoms.
6. Conclusion By adopting a modified version of EHT and by using the empirically determined parameters A, we have been able to obtain results for copper clusters which are broadly in agreement with ~~
(32) Tibbets, G. G.; Burkstrand, J. M.; Tracy, J. C. Phys. Reu. 8 1977, B15,3652. (33) Janak, J. F.; Williams, A. R.; Moruzzi, V.L. Phys. Reo. B 1975, Bll, 1522. (34) Bertoni, C. M.: Bisi, 0.;Calandra, C.; Manghi, F. Nuovo Cim. 1977, 388,96.
(35) Arlinghaus, F. J.; Gay, J. G.; Smith, J. R. Phys. Reu. B 1981, 823, 5152; 1980, B21,2201. (36) Citrin, P. H.; Wertheim, G. K.; Baer, V. Phys. Rev. Left.1978, 41, 1425. (37) Steiner, P.; Hufner, S.;Freeman, A . J.; Wang, D.Sulid Srure Commun. 1982, 44,619. (38) Egelhoff, W. F. Phys. Rev. 8 1984, 829,4769.
740 The Journal of Physical Chemistry, Vol. 92, No. 3, 1988
experiment. There are, however, some important exceptions and, in particular, the predicted conformations of small clusters are probably wrong and the EHT(A) estimates for the total energies of the larger clusters are much too small, showing no indication of convergence to the bulk cohesive energy. Following a referee’s suggestion, we find that much better values for the total energies of small and large clusters can be obtained using parameter set B although other properties are not well-predicted by these parameters. A possible explanation is that set A corresponds to a transition-state Hamiltonian particularly suitable for describing ionization and (but to a lesser extent) excitation processes whereas set B applies to an initial-state, Le., the ground-state, Hamiltonian. This is very speculative, however, and a simpler possibility for the difficulty we have found in getting satisfactory estimates of total energies and conformations is that the EHT energy formula (3) is a too crude and unreliable approximation. Because we use off-diagonal parameters which differ from the corresponding diagonal ones, we are able to obtain a d-band width which is much larger than that found in most EHT calculations and which, in fact, is similar to the bulk d-band width. It has been argued, on both experimental and theoretical grounds (see, for example, ref 39-44), that a real copper cluster consisting of, say, 50-100 atoms would have a rather narrow d band. If this were correct then, of course, our results would be wrong when applied to real clusters. However, the major reason for attempting cluster calculations of this sort is in order to investigate processes such as chemisorption on metal surfaces. For this to succeed it is necessary for the cluster to model the properties of the bulk or, at least, the properties of the metal surface. In either case, it would seem essential that the correct d-band width is reproduced. As we have shown, this can be done provided a modified version of EHT is used. Acknowledgment. We thank the Science and Engineering Research Council and ABM Chemicals Ltd., Leeds, for financial
Amos et al. support for S.A.M. Acknowledgments are due, also, to two referees for their interesting and helpful suggestions. Appendix If an orbital 4i is written as a linear combination of atomic orbitals w,
4i =
I
then the associated charge density will be P = Calaswlw, IS
assuming the coefficients and orbitals are real. On integrating over space we obtain
p = CarasSrs IS
where P , the total charge, will equal unity. We wish to divide P among the atomic orbitals, Le., to write p = CP, I
where pr is that part of P we associate with orbital w,. There is no unique choice for the pr, nor, indeed, a particularly obvious one because of the cross terms. Quite often Mulliken’s population analysis45is used and this divides the cross term arasSrsequally between orbitals r and s, so that PI
= a;
+ CarasSrs s#r
Unfortunately, there is no guarantee that the pr will be positive and, indeed, they often turn out to be negative. This can lead to physically meaningless results such as negative s contributions to an orbital. A way to avoid this difficulty is to ignore the cross terms by taking PI =
(39) Baetzold, R. C.; Mason, M. G.; Hamilton, J. F. J. Chem. Phys. 1980, 72, 366. (40) Roulet, H.; Mariot, J.-M.;Defour, G.; Hague, C. F. J . Phys. F 1980, 10, 1025. (41) Liang, K. S.; Salaneck, W. R.; Aksay, L. A. Solid State Commun. 1976, 19, 329. (42) Unwin, R.; Bradshaw, A. M . Chem. Phys. Lett. 1978, 58, 58. (43) Grunze, M. Chem. Phys. Lett. 1978, 58, 409. (44) Lee, S.-T.;Apai, G.; Mason, M. G.; Benbow, R.; Hurych, Z. Phys. Reu. E 1981, 823, 505.
Carw,
ar2/Car2
where the denominator has to be introduced to ensure the pr sum to unity. This is the formula we have used to evaluate s, p, and d DOS and surface and bulk DOS. Registry NO. CU,7440-50-8; CUI, 12190-70-4;C U ~66771-03-7; , Cud, 65357-62-2; C U ~66771-05-9; , C U ~681 , 12-39-0. (45) Mulliken, R. J . Chem. Phys. 1955, 23, 1833.
