J . Phys. Chem. 1989, 93, 1556-1565
1556
where all terms refer to 298 K and IE is the adiabatic ionization energy.& It is seen from Table VI that oxygen-containing bases always have higher ionization energies than their sulfur counterparts and AIE(R0,RS) varies from 1.10 to 1.44 eV. Moreover, the ionization energy differences are always greater than corresponding differences in the hydrogen-atom affinities, AHA, of oxide and sulfide ions. The AIE term therefore dominates the usual proton affinity order of oxide and sulfide bases and PA(RS) > PA(R0). Using our kinetically bracketed proton affinity of SiS, 170.3 f 2.0 kcal mol-’, and IE(SiS) = 10.26 eV7 we derive that the normal HA(SiS+) = 93.3 f 3.0 kcal mol-I is in good agreement with the average HA(RS+) = 88 f 5 kcal mol-l for the sulfurprotonated bases in Table VI. Moreover, AIE(Si0,SiS) = 1.35 eV7,45and is typical of the ionization energy differences among the sample of organic oxides and sulfides in Table VI. The heteroatom-site proton affinity order for S i 0 and SiS resolves, (46) Adiabatic ionization enthalpies for a sample of organic molecules at 350 K, IHpso,differ rarely from extensively compiled adiabatic ionization energies at 0 K, IEoo, by more than 0.05 eV.U Applying the stationary electron convention, it follows that the same difference should pertain to ionization energies at 350 and 0 K, and is likely to cancel in the term IEoo[(RO)-(RS)]. Therefore, IEoo[(RO)-(RS)] can substitute for IE0298[(RO)-(RS)] when only qualitative comparisons of the latter property are required or when approximate values HA[(RS+)-(RO+)] are derived from PA[(RS)-(RO)] .
therefore, to the exceptionally large HA(SiO+) = 143.4 f 2.8 kcal mol-’. This enthalpy is considerably larger than the average HA(RO+) = 108 f 12 kcal mol-’ or the empirical AH(CO+-H) = 112 kcal In fact, HA(Si0’) appears to be more representative of the metal oxide ions for which HA(Mg0’) = 133.7 f 18S;HA(CaO+) = 149.9 f 13.8, HA(Sr0’) = 140.7 f 9.2, and HA(BaO+) = 133.7 f 9.2 kcal mol-’ all at 0 K.43 The proton affinities of metal oxides range from 235 to 345 kcal mol-l I9 and the normal proton affinity of Si0 falls short of that range. Presumably the heteroatom-terminal proton affinity of silicon oxide reflects silicon’s intermediate position in the periodic table. Whether the relative order of heteroatom-site proton affinities for Si0 and SiS is representative of a general order for metal monoxides and monosulfides remains a topic for future experimental and theoretical investigations.
Acknowledgment. We thank the Natural Sciences and Engineering Research Council of Canada for financial support and Mr. John Houghton (Matheson Gases, Whitby, Ontario) who performed the purity analysis for allene. Registry No. H+, 12408-02-5; SiO, 10097-28-6; SiS, 113443-18-8; H O W , 66106-76-1; HSSi+, 76274-50-5; NH,, 7664-41-7; HCN, 7490-8; H2S, 7783-06-4; HzO,7732-18-5; C2H4,74-85-1; (CH3)20,11510-6; CH,COOH, 64-19-7; CH,CN, 75-05-8; CZHSOH, 64-17-5; HZCCCH2, 463-49-0; CH,OH, 67-56-1; HCOOH, 64-18-6; CO, 630-08-0; CS, 2944-05-0.
Corrected Effective Medium Method. 3. Application to Clusters of Mg and Cu Joel D. Kress,+ Mark S. Stave, and Andrew E. DePristo*,l Ames Laboratory-USDOE (Received: June 20, 1988)
and Department of Chemistry, Iowa State University, Ames, Iowa 5001 1
We present a number of new concepts in the recently developed corrected effective medium (CEM) theory. First, the general CEM theory is specialized to the case of an infinite periodic three-dimensional solid. Second, construction of a covalent embedding function of the electron density from diatomic and bulk binding potentials is detailed. On the basis of the smoothness of the function in interpolation between the (low density) diatomic and (high density) bulk limits, we argue that this function is dependent only on the types of atoms in a system and not on the number of atoms in the system (Le., a universal function of the electron density for particular types of atoms). Applications are made to the description of the structures and energies of homogeneous metal clusters containing from 2 to 201 atoms. Predictions of the equilibrium structures and energies of various high-symmetry structures of Mg,, Mg,, Cu,, Cu,, Cu,, CuI3, and C q 9 are compared to a number of ab initio and SCF-LSD results. The agreement is quite good even for the small clusters and becomes essentially perfect by Cu13. Selected two-dimensional distortions of four-, five-, and seven-atom clusters are considered in detail, with the energetics presented via contour plots. These enable the determination of barriers to isomerization and/or the existence of especially floppy modes of the metal cluster. Geometrical optimization is performed for various configurations of Mg,,, Cu,-,,, and CuI9. These accurate calculations on this range of clusters enable us to investigate the structural stability of the larger metal clusters for the first time. For still larger clusters of up to 201 atoms, we present binding energies for “spherical” bulk fragments to determine the onset of bulk binding.
I. Introduction Investigations of the structural and energetic characteristics of small metal and semiconductor clusters is an active area of theoretical’-20 and experimenta121-30research. Much work has been directed toward the calculation of the most stable geometries of small clusters in high-symmetry arrangements, AN with N i= 2-1 3, by reasonably accurate a b initio or first-principle methods. Larger clusters have been treated by cruder semiempirical or ab initio methods that are not capable of accurate energy calculations. In both cases, the very interesting question of the kinetic barriers to the isomerization reactions in such clusters has not attracted
’
Mobil Foundation and Henry Gilman predwtoral fellow. Present address: Group T-12,Mail Stop 5569, Los Alamos National Laboratory, Los Alamos, NM 87545. *Camilleand Henry Dreyfus Teacher-Scholar; Alfred P. Sloan fellow. 0022-3654/89/2093-1556$01.50/0
much attention; the degree of “floppiness” of various cluster modes has not been considered in detail; and studies of the many lower (1) See, for a general review; Koutecky, J.; Fantucci, P. Chem. Reu. 1986, 86, 539. (2) Pacchioni, G.; Koutecky, J. J . Chem. Phys. 1982, 77, 5850. (3) Jordan, K. D.; Simons J. J . Chem. Phys. 1980, 72, 2889. (4) Chiles, R. A.; Dykstra, C. E.; Jordan, K. D. J . Chem. Phys. 1981, 75,
1044. (5)
Delley, B.; Ellis, D. E.; Freeman, A.; Baerends, E. J.; Post, D. Phys. Reu. B 1983, 27, 2132. (6) Bachman, C.; Demuynck, J.; Veillard, A. Faraday Symp. R . Soc. 1W0, 14, 269. In Growth and Properties of Metal Clusters; Bourdon, J., Ed.; Elsevier: Amsterdam, 1980. (7) Flad, J.; Igel-Mann, G.; Preuss, H.; Stoll, H. Chem. Phys. 1984, 90,
257.
Jeung, G.H.; Pelissier, M.; Barthelat, J. C. Chem. Phys. Left. 1983, 369.
(8)
97,
0 1989 American Chemical Society
Corrected Effective Medium Method symmetry clusters have not appeared. There are two major reasons for this situation: (1) the accurate methods are too expensive computationally for the calculation of energies at many geometrical configurations; (2) the cruder methods are too inaccurate for such studies. In two recent we have derived and implemented a new approach to the calculation of interaction energies based upon explicit evaluation of the corrections to the effective medium theory33*34 using density functionals. The second referred to from now on as paper 11, presented a general formalism for the calculation of the interaction energy for any number and type of atoms (Ai,i = 1, ..., N). The acronym CEM-N is used to refer to this formalism. The basic idea of this theory is to replace each atom Ai interacting with all other atoms {A,,j = 1, ..., N, j # i) by atom Ai embedded in a spin-unpolarized jellium of density ni. The interaction or embedding energies for nearly all atoms through Cu are known as a function of the jellium density.35 The real system differs from the atom-jellium one due to inhomogeneities of the electron density and due to the point charge of the nuclei. We determine non-self-consistent corrections between the two systems, involving both Coulombic and kinetic-exchangecorrelation energies. The latter are minimized by proper identification of the homogeneous electron gas density corresponding to the real system. The details of the CEM-N theory along with illustrative tests for binding in a few first-row diatomic molecules (N2, C2, 0,) are contained in paper 11. (Results for H2, Liz, B2, C2, N2, 02, Rao, B. K.; Khanna, S. N.; Jena, P. Phys. Rev. E 1987, 36, 953. (IO) Ermler, W. C.; Kern, C. W.; Pitzer, R. M.; Winter, N. W. J . Chem. (9)
Phys. 1986,84, 3937.
( 1 1) Baykara, N. A.; Andzelm, J.; Salahub, D. R.; Baykara, S. Z. In?. J . Quantum Chem. 1986, 29, 1025. (12) Raghavachari, K. J . Chem. Phys. 1986,84, 5672. (13) Messmer, R. P. J . Vac. Sci. Technol. A 1984, 2, 899. (14) Upton, T. H.; Goddard, 111, W. A. Chemistry and Physics of Solid Surfaces: Vanselow, R., England, W., Eds.; CRC Press: Boca Raton, FL, 1985; Vol. 111. (15) Bauschlicher, Jr., C. W.; Pettersson, L. G. M. J . Chem. Phys. 1986,
84, 2226. (16) Pettersson,L. G. M.; Bauschlicher,Jr., C. W.; Halicioglu,T. J. Chem. Phys. 1987,87, 2205. (17) Upton, T. H. Phys. Reu. Lett. 1986,56,2168; J . Chem. Phys. 1987, 86, 7054. (18) Daudey, J. P.; Novaro, 0.;Kolos, W.; Berrondo, M. J. Chem. Phys. 1979, 71, 4297. (19) Post, D.; Baerends, E. J . Chem. Phys. Lett. 1982, 86, 176. (20) Wang, Y.; George, T. F.; Lindsay, D. M.; Beri, A. C. J . Chem. Phys. 1987,86, 3493. (21) See, for a number of general reviews: Morse, M. D. Chem. Reu. 1986, 86, 1049. Metal Clusters; Trager, F., zu Potlitz, G. Eds.; Springer: Berlin, 1986. The Physics and Chemistry of Small Clusters; Jena, P., Rao, B. K., Khanna, S. N., Eds.; Plenum: New York, 1987. (22) Morse, M. D.; Geusic, M. E.; Heath, J. R.; Smalley, R. E. J . Chem. Phys. 1985,83,2293. Morse, M. D.; Smalley, R. E. Eer. Bunsen-Ges Phys. Chem. 1984, 88, 228. (23) Hoffman 111, W. F.; Parks, E. K.; Nieman, G. C.; Pobo, L. G.; Riley,
The Journal of Physical Chemistry, Vol. 93, No. 4, 1989 1557 Fz, Naz, Mg2, A12, Si2, C12,and {K-Cu), have also been obtained in our lab that support those in paper 11.) For the present purposes, the important point is that the calculation of interaction energies with good accuracy at (relatively) low computational cost can be provided by this CEM-N theory. This sweeping statement must be qualified by the following fact: For near-quantitative accuracy, the SCF-LD embedding functions of Puska et al.35must be replaced by covalent embedding functions, which at this stage of the theoretical development, can be constructed only for homogeneous atomic clusters. The general problem of interpolation between the Puska and covalent embedding functions is especially severe for elements with large electronegativities. In addition, it is our opinion based upon a number of (unpublished) studies that in general metallic species will be treated more accurately than nonmetals. These qualifiers are meant to provide a measure of the current status of the C E M - N theory. Further developments may very well ease some of these problems, and ongoing work is directed toward this goal. However, there are a large number of interesting chemical systems that can be investigated by using the current level of the theory. The present paper provides examples in the structure and energetics of homogeneous metal clusters containing from 2 to 201 atoms. This paper is divided into four sections. A brief review of the CEM-N theory and its extension to infinite 3-D periodic systems is contained in section 11. Section I11 provides the results of CEM calculations on bulk and clusters of Mg and Cu, along with comparisons to other calculations. This section also contains a short discussion of the choice of embedding functions and detailed treatment of the extraction of the covalent embedding functions from diatomic and bulk binding potentials. A brief summary and conclusions are presented in section IV. 11. Theory In paper 11, we considered an N-body system of atoms ( A , i = 1, ..., NJ, where each A, can be any type of atom. We showed that the interaction energy, defined in obvious notation as the difference in energy between the interacting and noninteracting systems
AE({AjI)= E ( C A i ) - C E ( A J
(1)
could be rewritten exactly in the form
AE((A,))= ZAEJ(Ai;ni)+ AVc + AG({Ai/)
(2)
Equation 2 expresses the stabilization of the N-body cluster as a sum of three terms: (1) the sum of the embedding energies for each atom A , into jellium of density ni, AEJ(Ai;ni);(2) the difference in the Coulombic energy between the real system and all the atoms in jellium, AVc; (3) the difference in the sum of the kinetic, exchange, and correlation energies between the real system and all the atoms in jellium, AG({A,)).The latter term is given explicitly by
S . J. 2.Phvs. D. in Dress.
(24) Lii, K.; Pa&, E. K.; Richtsmeier, S. C.; Pobo, L. G.; Riley, S . J. J . Chem. Phys. 1985,83, 2882. (25) Ruatta, S. A,; Hanley, L.; Anderson, S. L. Chem. Phys. Lett. 1987, 129. 137: J . Chem. Phvs. 1987, 87. 260. (26) Kondow, T. J..Chem. Phys. 1987, 87, 3927. (27) Mandich, M. L.; Reents, Jr., W. D.; Bondybey, V. E. J . Phys. Chem. 1986, 90, 2315. Mandich, M. L.; Bondybey, V. E.; Reents, Jr., W. D. J . Chem. Phys. 1987,86, 4245. Reents, Jr., W. D.; Mujsce, A. M.; Bondybey, V. E.; Mandich, M. L. J . Chem. Phys. 1987, 86, 5568. (28) Whetten, R. L.; Cox, D. M.; Trevor, D. J.; Kaldor, A. Surf. Sci. 1985, 156, 8.
(29) O’Brien, S. C.; Liu, Y.; Zhang, Q.;Heath, J . R.; Tittel, F. K.; Curl, R. F.; Smalley, R. E. J . Chem. Phys. 1986, 84, 4074. (30) Jarrold, M. F.; Bower, J. E. J . Chem. Phys. 1986,85, 5373; 1987, 87, 1610; 1987, 87, 5728. Jarrold, M. F.; Bower, J. E.; Kraus, J. S. J . Chem.
Phys. 1987,86, 3876.
Kress, J. D.; DePristo, A. E. J . Chem. Phys. 1987, 87, 4700. Kress, J. D.; DePristo, A. E. J . Chem. Phys. 1988, 88, 2596. Norskov, J . K.; Lang, N. D. Phys. Reo. B 1980, 21, 2136. Stott, M. J.; Zaremba, E. Phys. Reu. B 1980, 22, 1564. (34) See: Jacobsen, K. N.; Norskov, J. K.; Puska, M. J. Phys. Reo. E 1987, 35, 7423, and references therein. (35) Puska, M. J.; Nieminen, R. M.; Manninen, M. Phys. Reo. E 1981, 24, 3037. Puska, M. J., private communication. (31) (32) (33)
where G ( A ) is the kinetic-exchange-correlation energy of the system A . The embedding energies are known quantities. In one case, they have been calculated for nearly all the atoms through Cu as a function of jellium density by using the SCF-LSD approach.35 These values are denoted as AEp(Ai;ni).In another case, they can be constructed from the binding potentials for the homonuclear diatomic and homogeneous bulk systems, as will be illustrated in section 111. These values are denoted as AEc(Ai;ni).In paper 11, we discussed the distinction between these two energy functions and showed that they are two special cases of a more general embedding function that is a function of both the jellium density and the work function of the jellium. Thus, the subscript C indicates that these embedding values are determined from essentially covalent binding interactions, in contrast to the SCF-LSD results, which have a large contribution from ionic binding. These covalent embedding values provide interaction energies of chemical accuracy that is not possible by using the Puska et al. embedding
1558 The Journal of Physical Chemistry, Vol. 93, No. 4, 1989 function. This point will also be illustrated in detail in section 111. However, the qualitative behavior of binding energy variation with internuclear separation is often similar by using the Puska et al. and covalent embedding functions. Evaluating eq 2 via S C F techniques would be as difficult as evaluating the original eq 1. Thus, in the CEM approach we utilize the approximation of superposition of electron densities, or that the electron density at any point in space, r, is the sum of the densities from each atom: n(r) = Zn(r;Ai)
(4)
This superposition of atomic densities is the major approximation in implementation of the theory, but it is utilized only to calculate the correction terms: AVc and AG, Since this non-self-consistent part of the C E M - N formalism is not expected to be as accurate as the self-consistent part, we minimize AG with respect to the choices of the jellium densities Ini). This yields the general expression for the jellium density, in terms of the atomic up- and down-spin densities (n+(Ai), n-(Ai)):
ni =
c(1 - 6ij)${ni(Ai)n'(Aj) + n-(Ai)n-(Aj)) d r / Z , J
(5a)
which reduces for unpolarized densities to
where Zi is the atomic number. The superposition approximation makes the Coulombic difference simply the sum of the atom-atom Coulomb interactions. Finally, we note the functionals used in evaluation of AG. The kinetic energy functional is an accurate Pade a p p r ~ x i m a t i o nin~ ~ IVnl/n4/3, which approximately sums the full series in the gradients of the electron density. The local Dirac3' exchange functional and the local Gunnarsson-Lundqvist (GL) correlation functional3* are used. This finishes the brief outline of the CEM-N formalism. The reader interested in more details should consult papers I and 11. However, it is necessary to understand only the basic ideas discussed above to follow the derivation of the bulk limit of these equations. We now proceed to derive the N m limit in eq 2 and provide the C E M expression for the cohesive energy for an infinite solid with 3-D translational periodicity. Under this condition, each Ai is in an identical environment. Thus, we can focus on a single atom Ai in the lattice and reexpress eq 2 as
-
-
where Vc(ij) is the Coulomb interaction between atoms i and j . The limit of eq 6 as N must be taken to determine the cohesive energy (Le., the interaction energy per atom). This can be accomplished for the last term in eq 6 most readily in terms of a quadrature scheme introduced in the appendix of paper 11: reexpressing a multicenter integral of one global function by a multicenter integral over many local functions, defined by the distance to each atomic center. This allows for a reexpression of the multicenter integral with respect to integrals around each center, yielding where the WS(i) indicates that the integration of AG extends only over thefinite volume of a Wigner-Seitz cell inscribed about the atom Ai. Substituting eq 7 into eq 6 yields the CEM-N cohesive energy expression AE(IAiJ)/N = AEJ(Ai?i)
+ Y*X Vc(i4) + AG(IAj))ws(il j#r
The indexj in eq 5-7 effectively sums over an infinite number of atoms. In practice, the summation in eq 5 converges very quickly because n(r;A,) decays exponentially in r = Irl, and thus any atom A, situated at a distance away from Ai where n(r;A,) effectively equals zero can be neglected in the sum. The Coulomb integral extends to a longer range, but the summation does not cause any problem for neutral atoms. If the lattice contained ions, then special care would have to be used in converging this summation. Next, we discuss a couple of computational details of the CEM approach. We have constructed an even-tempered Gaussian basis39 to represent the atomic densities n(r;Ai) which are generated from Slater-type atomic Hartree-Fock densities.40 Use of Gaussians allows for efficient analytic evaluation41 of the Vc(ij) and the density overlaps that are combined to form the jellium densities in eq 5. Thus, the AG term in eq 2 and 8 is the only term that is calculated by using a numerical quadrature scheme. A radial-angular grid is centered about each Ai, and AG is evaluated by using the quadrature scheme discussed earlier and detailed in paper 11. Specifically (1) the radial integration is performed by using an extended Bode's rule42( 2 ) the 0 integration is performed by using a Gauss-Legendre and (3) the 4 integration is performed by using a Gauss-Chebyshev q ~ a d r a t u r e . We ~~ should emphasize that more optimal integration schemes43can be implemented in this part of the calculation as they become available. Since the major computational time occurs in this numerical integration scheme, a more efficient integration scheme will lower the computational time proportionally. 111. Results Diatomics, clusters, and bulk structures of both Mg and Cu are considered in this paper. We must determine the spatial symmetry and spin polarization of the atomic densities that are superposed to form the cluster and/or bulk density (both topics of which were investigated thoroughly in paper I1 for the first-row homonuclear diatomic molecules). Here, this poses no problems since Mg is a closed-shell atom, which implies a spherically symmetric and spin-unpolarized atomic density. For the Cu atom we used the Hartree-Fock valence shell that is nearly closed with a configuration of (3d9,4s2) as tabulated in ref 40. As pointed out in paper I, this is clearly not the valence configuration of the bulk Cu band, which is known to have a partially filled s band, suggesting a valence configuration of (3dI0,4s'). Yet, as we pointed out in paper 11, the choice of valence configuration used within the superposition approximation is a minor consequence, since we are only trying to model the gross features of the many-body system. The decisive factor in favor of the (3d*,4s2) over the (3d"+',4si) configuration for all the transition metals we have considered to date (Sc through Cu) concerns the construction of the covalent semiempirical embedding functions. In general, we have found the two diffuse 4s electrons of the former allows us to match the diatomic and bulk contributions to the embedding functions more smoothly. This will become clear in part B. As a further approximation, we have set the 3d shell spherically symmetric and spin-unpolarized for the Cu atoms. This has a negligible effect. A . Bulk Cohesive Energies. For Mg and Cu, we have calculated the cohesive potentials using the self-consistent Puska et al. embedding energies for AEJ in eq 8. The results for the various lattice structures considered for each metal are presented in Table I, along with the experimentally observed lattice structures and values. The hcp calculation was performed by constraining the ratio of the c / u lattice constants to the ratio that is experimentally observed,44with the potential then predicted as a function of u.
(8)
The jellium density is constructed as for the finite cluster in eq 5 but now need only be done once because all atoms are identical. (36) DePristo, A . E.; Kress, J. D. Phys. Rev. A 1987, 35, 438. (37) Dirac, P.A. M. Proc. Cambridge Philos. SOC.1930, 26, 316. (38) Gunnarsson, 0.; Lundqvist, B. I . Phys. Reo. B 1976, 13, 4274.
Kress et al.
(39) Schmidt, M. W.; Ruedenberg, K. J . Chem. Phys. 1979, 71, 3951. (40) Clementi, E. I B M J . Res. Deu. Suppl. 1965, 9. Bagus, P. S . ; Gilbert, T. L.; Roothan, C. J. J . Chem. Phys. 1972, 56, 5195. (41) Huzinaga, S. Prog. Theor. Phys. Suppl. 1967, 40, 279. (42) Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Func(ions; Dover: New York, 1972. (43) (a) Boerrigter, P. M.; te Velde, G.; Baerends, E. J. fnr. J . Quanl. Chem. 1988 33, 87. (b) Becke, A. D. J . Chem. Phys. 1988, 88, 2547.
The Journal of Physical Chemistry, VOI.93, No. 4, 1989 1559
Corrected Effective Medium Method TABLE I: CEM-N Cohesive Potentials for Various Lattices of Mg and Cu Metals Calculated with the Puska et at. Embedding Functiod5 lattice lattice NND,O cohesive metal type const, bohrs bohrs energy, eV Mg hcpb 5.76 5.76 1.62 fcc 8.16 5.77 1.58 bcc 6.44 5.59 1.57 (hCPY (6.07) (6.07) (1.53) cu fcc 6.72 4.75 4.26 bcc 5.31 4.65 4.30 (fee) (6.82) (4.82) (3.50)
"NND is the nearest-neighbor distance for the given lattice type and constant. bThe potential for the hcp lattice was determined by fixing the ratio of the c / a lattice constants to that observed experimentally.44 The a lattice constant is given. C I n parentheses are the experimentally observed parameters as found in ref 44.
Examining Mg first in Table I, we note that all three lattice structures are predicted to be of essentially equal stability with only slight preference for the hcp value. For the hcp lattice, the cohesive energy is overestimated by 5.9% while the lattice constant is underestimated by 5.1% with respect to experiment." Examining Cu next, we note that both lattice structures are predicted again to be of essentially equal stability, but now for the fcc lattice the cohesive energy is overestimated by 22% while the lattice constant is underestimated by 1.5% with respect to experiment." This is a general characteristic that we have found for essentially all bulk systems including semiconductors and insulator^.^^^^^ These comparisons demonstrate that the use of the Puska et al. embedding functions provides a reasonable qualitative description of the cohesive energies and perhaps even semiquantitative geometrical spacings for these metals within the superposition of atomic densities approximation as applied within the CEM-N formalism. An alternative embedding function must be considered to provide quantitative binding energies. The central contribution this paper makes toward this goal is to put forth the necessary ingredients to construct a semiempirical embedding function from the knowledge of the experimental binding potentials for only the homonuclear diatomic molecules ( N = 2) and bulk solids (N = a). With the two extremes of N-atom clusters described accurately by definition, we then propose that the new semiempirical embedding function is universal with respect to N for a given element as applied within the C E M approach. This leaves for our consideration N-atom systems both small (Le., metal clusters) and large (Le., semiinfinite metal surfaces). Such properties as equilibrium geometries and interaction energies for the former systems will be described in this paper, whereas the treatment of metal surfaces is deferred to a later paper.46 B . Covalent Embedding Energy Functions. Semiempirical embedding functions were presented in paper I1 for atoms in the first row of the periodic table. These were constructed solely by inverting the CEM-N binding potentials with respect to the experimental homonuclear diatomic binding curves.47 These embedding functions were designated as covalent since they were derived from homonuclear diatomics, which are assumed to bind covalently. Specifically, these covalent embedding functions AE,(A;n) were derived from eq 2 as
AEc(A;n)= (AE(A2)- AVc - A G ( A z ) ) / 2
(9)
with the Morse potential for A 2 substituted into the right-hand side of eq 9. A new type of covalent embedding function, proposed earlier in this paper, is constructed by also inverting the CEM-N cohesive potential for bulk solids with respect to experiment and smoothly combining this information with the information derived from eq 9. In our analysis, the bulk experimental potential is constructed ~
(44) Kittel, C. Introduction to Solid Stare Physics, 4th ed.; Wiley: New York, 1971. (45) Raeker, T.; DePristo, A. E. Unpublished results. (46) Raeker, T ; DePristo, A. E. To be submitted to Phys. Reu. B. (47) Huber, H. P.; Herzberg, G. Constants of Diatomic Molecules; Van Nostrand. New York, 1979
n
x 10-3
n
(au)
b
v m
>
2
x
'0-3 (CUI
l
INTERNUCLERR DISTANCE (au)
BONO LENGTH (au)
Figure 1. Experimental potentials for Mg from which the Mg semiempirical covalent embedding function is derived, plotted as a function
of both the internuclear spacing (nearest-neighbordistance for the bulk; bond length for the diatomic) and the corresponding sampled jellium density n: (a) bulk metal; (b) diatomic molecule. MS Emaecid i -3 = u i c I TS
6 -
se f 3orsis:er:
3
4
*
CLUSI-' 3 1 L. Figure 2. Semiempirical covalent embedding function for Mg. Also shown is the SCF-LD embedding function as calculated by Puska et aLss The embedding energy is plotted as a function of jellium density. E_tCTFI
I
4
by considering a harmonic expansion about the equilibrium lattice constant, with the bulk modulus (3.54 X 10" dyn/cm2 for Mg and 1.37 X lo1*dyn/cm2 for Cu)" providing the second derivative of the potential. The contribution to the covalent embedding function from the bulk is then derived by inverting eq 8 upon substitution of the experimental potential into the left-hand side. To illustrate this procedure in detail, we provide a step-by-step analysis for the construction of the covalent embedding energy for Mg. In Figure l a , we plot the harmonic cohesive potential from explicit points at 90%, 95%, loo%, 105%, and 110% of the equilibrium lattice spacing. In Figure lb, we show the Morse potential for the diatomic. The former will be less accurate at large expansions while the latter will be less accurate at small bond lengths. Note that the potentials are plotted as functions of both the internuclear distance and the corresponding jellium density n. This information is inverted as prescribed to obtain AEc(Mg;n) shown in Figure 2. The four highest density points of hEc(Mg;n) in Figure 2 correspond to the four highest density points in Figure la. To interpolate smoothly between the high- and low-density limits, we truncated the bulk contribution to AEc(Mg;n) at the 105% point denoted by the arrow on Figure la. The point (filled square) denoted "N = 2 last h o t " on Figure 2 corresponds to the point on the curve in Figure lb, with the still lower density points (not shown by explicit symbols) corresponding to the portion of the potential in Figure 1b starting at the (filled square) point and extending out to large values of bond length. Upon examination of both Figures 1 and 2, one notes right away the difference in magnitude between the values of n derived from the bulk vs diatomic systems. This is due to the fact that a bulk Mg atom samples (mainly) six nearest neighbors and six (very close to first but) second nearest neighbors to provide a jellium density at a given internuclear separation, whereas the Mg atom in the diatomic samples only one atom. The significant variation in magnitudes provides a sensitive test of the universality of the embedding function via the degree of smoothness in interpolation
The Journal of Physical Chemistry, Vol. 93, No. 4, 1989
-
x '3.3
(au;
-
1.3 0.50 0.12 0.03
15.6
11.4
3.2
6.0
1,
,
'
'
INTERNUCLEAR DISTANCE
4.3
10-3
,I
BONO -ENSTH
(OL 1
(OL:
Figure 3. Experimental potentials for Cu from which the Cu semiempirical covalent embedding function is derived, plotted as a function of both the internuclear spacing (nearest-neighbordistance for the bulk; bond length for the diatomic) and the corresponding sampled jellium density n: (a) bulk metal; (b) diatomic molecule.
-z-r-TFi -
__
,
I X E -I I3 Z Figure 4. Semiempirical covalent embedding function for Cu. Also shown is the SCF-LD embedding function as calculated by Puska et al.35 The embedding energy is plotted as a function of jellium density. \
between the two different regions. Clearly, A&(Mg;n) is a smooth function. Note also that AEc(Mg;n) is in rather good general agreement with AEp(Mg;n), providing essentially a smoothing of the slight waviness of the Puska et al. SCF-LD values. Of course, note also that the Puska et al. calculations do not extend to low enough density to determine the diatomic embedding region. In a manner analogous to the Mg case, we have constructed the covalent embedding function for Cu. In Figure 3a , we plot the experimental harmonic cohesive potential for the metal as a function of the nearest-neighbor distance (NND) and as a function of jellium density. Likewise, in Figure 3b we plot the Morse potential for the diatomic molecule as a function of bond length and of jellium density. Note again the discrepancy in size of the jellium densities derived from the bulk vs the diatomic, although the difference is not as large as in the Mg case. We were forced to truncate the contribution from the diatomic potential closer to the equilibrium bond length [see the (filled square) point in Figure 3b] than in the Mg case to smoothly match the two contributions. Additionally, the bulk contribiution was again truncated at the point denoted by the arrow on Figure 3a. The net result is the covalent embedding function for Cu as shown in Figure 4. In contrast to the Mg case, the SCF-LD Puska et al. embedding function (solid line on Figure 4) is below the covalent curve in the region of the four knots provided by the bulk potential. This behavior is analogous to that found for homonuclear diatomic molecules in paper I1 and can be traced to the overestimation of ionic bonding character in homonuclear species when using the SCF-LD atom-jellium embedding energies at high jellium densities. A detailed discussion of this point can also be found in ref 46. The elimination of the 110% bulk point for both Mg and Cu is in accord with the increased stiffness of their bulk potentials. Analogously, the short-range repulsive Morse potential is likely to be too repulsive for a transition-metal dimer due to interchange of 4s and 3d electrons. The differences are quite small between
Kress et al. the CEM-N predicted Cu2 potential and the Morse curve in any case. For example, the use of the AEc(Cu;n) leads to a hardsphere radius of 3.09 bohrs while the Morse potential yields 3.26 bohrs, and the two potentials parallel each other in the region of the repulsive wall. It is perhaps worthwhile to note the insight provided by the C E M - N theory into the bonding relationship between the dimer and bulk system. Both utilize the same embedding function, albeit at rather different densities, with the increased relative stability of the bulk system being due to a larger Coulomb energy and a lower correction energy. The latter is a consequence of the more homogeneous electron density. C . Cluster Geometries. Before considering the predicted CEM-N potentials for small clusters of Mg and Cu atoms, we must describe (1) the various geometries that were considered for each value of N and (2) the coordinates that were varied. When possible, we utilize the notation in ref 1: [n.m] standing for n atoms in the mth geometrical arrangement given in chart 1 of ref 1. We also specify the point group when it is useful. Line drawings of the various structures are shown when a verbal description is not clear and/or when these structures are not shown in ref 1. The structures are described as follows: 3 atoms, linear (D-h): two bonds of length a. 3 atoms, bent ([3.2], C2"):two edges of a with an angle @ between them. 4 atoms, linear (Om,,): inner bond of a and outer bond of b. 4 atoms, T-shaped ([4.5], C2u):bonds of a and b shown in Figure 5. 4 atoms, rhombus ([4.1], D2h): all edges of a with an acute angle 0. 4 atoms, tetrahedron ([4.3], Td): N N D of a. 5 atoms, linear ( D m h ) : two inner bonds of length a and two outer bonds of length b. 5 atoms, centered tetrahedron ([5.6], T d ) : N N D of a. 5 atoms, square pyramid ([5.3], C&): square's edges of a and length from apex atom to square's atoms of b. 5 atoms, trigonal bipyramid (15.21, &): triangle's sides of a and length from apex atoms to triangle's atoms of b. 6 atoms, dicapped tetrahedron ([6.5], Cb): N N D of a shown in Figure 5. 6 atoms, octahedron ([6.3], oh): N N D of a. 7 atoms, tricapped tetrahedron ([7.2], C3J: N N D of a shown in Figure 5. 7 atoms, fccl: additional atom at corner 1 of 6 atoms at the faces of a cube as shown in Figure 5. 7 atoms, pentagonal bipyramid ([7.1], DSh):pentagon's sides of a and length from apex atoms to pentagon's atoms of b. 8 atoms, cubic ([8.3], Oh):N N D of a. 8 atoms, quadcapped tetrahedron ([8.1], Td): N N D of a shown in Figure 5. 8 atoms, square antiprism (CdU):N N D of a in Figure 5. 8 atoms, fccl2: additional atoms at corners 1 and 2 in Figure 5. 8 atoms, fccl3: additional atoms at corners 1 and 3 in Figure 5. 8 atoms, fccl7: additidnal atoms at corners 1 and 7 in Figure 5. 9 atoms, fcc123: additional atoms at corners 1, 2, and 3 in Figure 5. 9 atoms, fcc127: additional atoms at corners 1, 2, and 7 in Figure 5. 9 atoms, fcc136: additional atoms at corners 1, 3, and 6 in Figure 5. 10 atoms, fcc1234: additional atoms at corners 1, 2, 3, and 4 in Figure 5. 10 atoms, fcc1235: additional atoms a t corners 1, 2, 3, and 5 in Figure 5. 10 atoms, fcc1236: additional atoms at corners 1, 2, 3, and 6 in Figure 5. 10 atoms, fcc1238: additional atoms at corners 1, 2, 3, and 8 in Figure 5. 10 atoms, fcc1357: additional atoms at corners 1, 3, 5, and 7 in Figure 5.
The Journal of Physical Chemistry, Vol. 93, No. 4, 1989 1561
Corrected Effective Medium Method
TABLE 11: Binding Characteristics for Mg Clusters Calculated with the Covalent Embedding Function a,'
C2v 6
c2v 4
c4v
8
Ih 13
FCC12
FCC123
FCC127
bohrs 7.35 5.33 3 linear 5.77 3 bent 5.82 5.22 4 linear 5.72 4 rhombus 5.81 4 tetrahedron 5.88 4.89 6.4 6.4 7.75 6.3 5 linear 5.71 5 centered tetrahedron 5.39 5 sq pyram 5.78 5 trigonal bipyram 5.92 6 dicapped tetrahedron 5.85 6 octahedron 5.71 m bulk (hcp)C 6.07
Td 8
D3h 13
FCCl
N geometry 2 lineard
c3v 7
Oh 13
FCCl3
5.77 bohrs 5.91 bohrs 5.85 bohrs 9.85 bohrs
method CEM-N SCF-CI CEM-N CEM-N SCF-CI CEM-N CEM-N CEM-N ~ 0 . 6 8 SCF-CI 0.18 SCF-CI 0.11 ECP-CI 0.03 SCF 0.26 CEPA 0.16 CEM-N 0.21 CEM-N 0.39 CEM-N 0.40 CEM-N 0.47 CEM-N 0.49 CEM-N 1.51 CEM-N
ref 9 9
9 4 2 3 4
"The a parameter for most geometries is defined in Figure 5 or Chart I of ref 1. bIf applicable, either the bond distance, b, or the bond angle, @, is defined in Figure 5. BEIN is the binding energy per atom. dThe experimental potential is predicted by definition, with the binding parameters obtained from ref 47. CTheexperimental potential is predicted by definition, with the binding parameters obtained from ref 44. The a parameter is the a lattice constant (which is also the NND), and the b parameter is the c lattice constant.
FCC17
FCC136
BEIN: eV 0.027 ~0.45 0.09 59.7O 0.18 60.0° ~0.5 5.77 bohrs 0.13 62.5' 0.27 0.32 b or pb
13 atoms, icosahedron ([13.1], Zh): N N D of a shown in Figure
FCC1234
5.
FCC1235
FCC45678
FCC245678
FCC1236
FCC1238
FCC34568
FCC234568
FCC1357
FCC24578
FCC2345678
FCC1368
FsC345678
-
FCC12345678
Figure 5. Definition of the geometric binding parameters for some of the
clusters. 10 atoms, fcc1368: additional atoms at corners 1, 3, 6, and 8 in Figure 5. 11 atoms, fcc45678: additional atoms at corners 4, 5, 6,7, and 8 in Figure 5. 11 atoms, fcc34568: additional atoms at corners 3,4, 5, 6, and 8 in Figure 5. 11 atoms, fcc24578: additional atoms at corners 2,4, 5, 7, and 8 in Figure 5. 12 atoms, fcc345678: additional atoms at corners 3, 4, 5, 6, 7, and 8 in Figure 5. 12 atoms, fcc245678: additional atoms at corners 2, 4, 5 , 6, 7, and 8 in Figure 5. 12 atoms, fcc234568: additional atoms at corners 2, 3, 4, 5, 6, and 8 in Figure 5. 13 atoms, fcc2345678: additional atoms at corners 2, 3, 4, 5, 6, 7, and 8 in Figure 5. 13 atoms, hcp structure ( [ 13.31, D3h): N N D of a shown in Figure 5. 13 atoms, cubooctahedron ([ 13.21, oh):N N D of a shown in Figure 5.
14 atoms, fcc12345678: additional atoms at corners 1, 2, 3, 4, 5, 6, 7, and 8 in Figure 5. 19 atoms, capped cubooctahedron ( o h ) : N N D of a . D. M g Clusters. In this section we consider the CEM-N potentials for Mg clusters generated by using the semiempirical covalent embedding function. By definition the potential for the dimer and the cohesive energy of the metal are reproduced exactly within the CEM-N approach. For convenience we provide the experimental binding energy per atom (BEIN) and equilibrium bond length for N = 2 and the experimental cohesive energy and N N D for N = a in Table 11, along with all the binding energies and distances for the various clusters. The number of integration points was sufficient to ensure a precision of i= 0.01 eV in energies and i= 0.03 bohr in equilibrium distances. These results are discussed below. The calculations for a linear three-atom cluster yielded a B E I N over 3 times larger than that for the dimer, and a very large contraction of the internuclear distance to 5.77 bohrs from 7.35 bohrs for the dimer. When p was varied also ( N = 3, bent in Table 11), the equilibrium arrangement had a = 5.82 bohrs and p = 59.7", almost a perfect equilateral triangle. Such a structure maximizes the number of nearest neighbors for a three-atom system, with the additional nearest neighbor providing for twice the B E I N with respect to the linear configuration. Interestingly, the bond lengths are about the same for the linear and bent configurations. The doubling of B E I N between the linear (two N N pairs) and triangular (three N N pairs) configurations 'indicates that a pairwise additive potential model would be rather poor. To examine the four-atom system in a logical manner, we first considered the linear arrangement where the CEM-N values were a = 5.72 bohrs, b = 5.77 bohrs, and B E I N = 0.13 eV, as listed in Table 11. Among the examined planar systems, the rhombus was the most stable with a = 5.81 bohrs, /3 = 62.5O, and B E I N = 0.27 eV. This geometry is more than twice as stable as the linear structure, again due to the increased number of nearest neighbors provided by the rhombus. Next allowing deformation to a three-dimensional structure, the tetrahedron, we obtained N N D = 5.88 bohrs, and B E I N = 0.32 eV, a gain of 0.05 eV in BEIN
Kress et al.
1562 The Journal of Physical Chemistry, Vol. 93, No. 4, 1989
t
Mg4: RHOMBUS TO TETRAHEDRAL
70.0
60.0
50.0 I\
a
:
"
40.0
*
LL
30.0
20.0
10.0
L---L-_iI-i---L-
~
i _
NND (au)
Figure 6. Contour plot of BE for the symmetricaldistortion of M& from a rhombus (4 = Oo) to a tetrahedron (4 = 57.74').
and an expansion of the NND. A more detailed investigation of the deformation is provided by the contour plot in Figure 6, in which the two atoms on the long diagonal of the rhombus were allowed to distort out of the plane by an angle q5 (=54.74' for a tetrahedron). There is clearly no barrier to this deformation even though the rhombus and tetrahedral structures are reasonably close in BE. The increase in BE occurs rapidly rather late in the deformation. The tetrahedral geometry is the most stable since it provides the largest possible number of nearest neighbor interactions in a four atom system. For the five-atom Mg system, the linear arrangement was calculated to yield a = 5.71 bohrs, b = 5.77 bohrs, and B E I N = 0.16 eV. This increase in B E I N with respect to the four-atom linear cluster is consistent with the increase observed from the three- to four-atom linear cluster. Relying on experience with the four-atom clusters, only three-dimensional structures were examined with the goal of determining the most stable structure for a given value of N. The equilibrium for the centered tetrahedron had a = 5.39 bohr and B E I N = 0.21 eV. This is only slightly more stable than the linear structure since there is no special stability associated with 4-fold coordination of a central Mg atom. In contrast to the case for the Td structure of the four-atom system, it is not possible to construct a five-atom structure where each atom has the same coordination, but two structures come close: the square pyramid and the trigonal bipyramid. For the square pyramid we determined a = 5.78 bohrs, b = 5.91 bohrs, and B E I N = 0.39 eV. The trigonal bipyramid was slightly more stable, with a = 5.92 bohrs, 6 = 5.85 bohrs, and B E I N = 0.40 eV. This slight increase in stability is due to the fact that the trigonal bipyramid provides a greater average coordination than the square pyramidal structure. This most stable structure was also investigated for symmetrical distortions as shown in Figure 7, where the abscissa is a and the ordinate is b. The contours indicate that the molecule is quite floppy, with the contour at BE = -1.95 eV allowing for a change of -1 bohr in both bonds! This subject will be addressed further in the cases of the D3h five-atom and the DShseven-atom Cu clusters. For the six-atom cluster, we considered two high-symmetry three-dimensional structures: the dicapped tetrahedron and the octahedron. Note the similarity in BEIN. The deformation of one structure into another was not investigated since a minimum energy path could not be determined in any simple manner. This would be extremely interesting to study by using a b initio or first-principle methods with analytic first and second derivatives. At this point let us summarize the geometry-optimized CEM-N potentials presented thus far for Mg. The most stable arrangements found for each respective value of N were (nearly) equilateral triangle (3), tetrahedron (4), trigonal bipyramid (9,and octahedron (6). Generally, C E M - N has predicted the geometry
Figure 7. Contour plot of BE for the variation of both the axial-toequatorial and the equatorial-to-equatorialdistances in Mg5. The arrows indicate the value of the axial-to-axial distance at the selected points
with the largest average coordination as the most stable for these clusters. By referring to Table 11, we note the following trends in the potentials: The internuclear atomic distances are all contracted significantly compared to the dimer bond distance and are all somewhat smaller than the bulk NND. Furthermore, the dominant internuclear distance decreases slightly as N is increased from 5 to 6. At some value of N , the trend must reverse to m limit. By contrast reproduce the bulk structure in the N B E I N increases uniformly as the cluster size increases, although it is still much smaller than the bulk value at N = 6. Although many HF-SCF c a l c u l a t i ~ n have s ~ ~ been ~ ~ ~performed ~ for small Mg clusters, they have predicted either an unbound or a very weakly bound cluster with respect to the free Mg atoms. Only by increasing the effort of the calculation and providing for electron correlation is a bound cluster found. Specifically, the potential for a tetrahedral cluster of Mg atoms has been generated by using a second-order coupled electron pair approximation (CEPA-2).4 This calculation yielded an equilibrium N N D = 6.3 bohrs and BEIN = 0.26 eV. Although no experimental data exist for Mg, clusters, the agreement between the CEPA-2 and the CEM-N potentials is reasonable. Additionally, a self-consistent electron pairs method (SCEP) yielded a N N D = 6.4. The discrepancy between these calculations and the CEM-N N N D of 5.88 bohrs reflects the difficulty in assigning an equilibrium N N D to a weakly bound cluster such as Mg4. In comparison to the B E I N predicted by the SCEP method4 (=0.177 eV) and by an effective core potential-multireference CI approach with Davidson's correction* (=O.1 13 eV), the C E M - N value for the B E I N is very good, especially since the above approaches tend to underestimate the binding for other metallic cluster systems as we shall see for Cu. Note that the SCF-LCAO-MO-MP ( M P = Moeller Plesset) results from ref 9 yield uniformly small bond lengths and large BEIN, and must be considered unreliable due to the poor predictions for the dimer.48 Clusters of Mg atoms containing more than six atoms were also studied, although geometry optimizations were not performed. Instead, the potentials were calculated with geometries consisting of "spherically" symmetric fragments of the bulk metal lattice with atomic distances fixed to be those observed experimentally in the bulk. By spherical, we mean that we inscribed a given radius about an atom in the periodic lattice and included all other atoms contained in the sphere. By increasing the size of the sphere, we can increase the size of the corresponding spherical cluster. The B E I N for such Mg clusters of up to 201 atoms are presented in Figure 8 (along with the geometrically optimized values for N = 2-6) as a function of X = 1 - (1 - N-'/3)3, as suggested by
-
(48) This inaccuracy is a reflection of the limited basis set size used in the calculations of ref 9 according to P. Jena, private communication.
Corrected Effective Medium Method
MgN v
The Journal of Physical Chemistry, Vol. 93, No. 4, 1989 1563
C I us t e r s
"1 rLK
TABLE 111: Binding Characteristics for Cu Clusters Calculated with the Covalent Embedding Function a," BEIN: N geometry bohrs b or fib eV method ref 2 lineard 1.02 CEM-N 4.20 0.94 ECP-LC 7 4.30 3 linear 1.14 CEM-N 4.34 0.86 ECP-LC 7 4.44 0.51 ECP-CI 8 4.82 0.43 SCF 6 4.44 3 planar 1.32 CEM-N 4.44 62.2'
1.5
cr &
201H0
'.2
135
c
i A
c o 3 C
0 2
x
=
C 6
0 4 .
I
.N
3.6
.G
1/31:
Figure 8. CEM-N binding energy per atom for N-atom Mg clusters constructed as fixed "spherical" fragments of the bulk hcp lattice (see text) for N 1 13 and geometrically optimized for N = 2-6. All calculations used the covalent embedding function. The values are plotted as a function of the fractional surface volume, X.
Stwalley and de Llano.49 This variable provides an estimate of the fractional surface volume (1 monolayer thick) of a cluster containing N atoms. It has been suggested" that X is the variable one should use when attempting (linear) extrapolation of data on finite clusters to that of the bulk system. A test of this hypothesis is provided in Figure 8 in an attempt to predict the CEM-N bulk cohesive energy from the knowledge of the CEM-N binding energies per atom for finite clusters of M g atoms. Assuming that fixing the internuclear spacings of the bulk lattice is a good approximation to the corresponding distances in the clusters for N 1 13, we note right away that a linear extrapolation in X provides a bulk BEIN 1.9 eV, which overestimates the actual value by about -25%. If the smallest clusters were not geometrically optimized, the curve would be even less linear. E . Cu Clusters. Next, we consider the CEM-N potentials for Cu clusters generated using the semiempirical covalent embedding function. Cu clusters are advantageous to study since there are a number of equilibrium binding energies and geometries that exist' to which we can compare the CEM-N results. Because of the accuracy of the CEM-N theory for metals, we compare directly only to high quality ab initio or first-principle calculations,'-8 but the interested reader may consult Table VI1 of ref 1 for a multitude of other results based upon diatomics-in-molecules (DIM), complete neglect of differential overlap (CNDO), and extended Huckel (EH) techniques. Again, the diatomic and bulk potentials are reproduced exactly by definition with the CEM-N approach. The experimental binding characteristics for both are provided in Table I11 for reference. The first system considered is the linear three-atom structure, for which we found a = 4.34 bohrs and BEIN = 1.14 eV. This is only slightly more stable than the diatomic BEIN, which should be contrasted with the 3-fold increase seen between Mg, and Mg,. When 0 is allowed to vary, allowing for a planar structure, the predictions are a = 4.44 bohrs, /3 = 62.2O, and BEIN = 1.32 eV. This geometry forms almost an equilateral triangle, and as in the M g case, this triangular structure is more stable than the linear structure due to the increased number of nearest neighbors. But here the triangular structure is only slightly more stable than the linear rather than nearly twice as stable as for the Mg, system. In comparison to potentials calculated with other methods, an effective core potential with a local correlation energy correction (ECP-LC) approach' yields BEIN = 0.89 eV and a discrete variable (DV)-Xa c a l c ~ l a t i o nyields '~ BEIN = 1.06 eV (for an unoptimized a = 4.82 bohrs the bulk NND), with both studies considering a triangular structure. The CEM-N predicted BEIN is larger than either of the above predictions. For the four-atom system, we first considered a linear structure which yielded a = 4.47 bohrs, b = 4.33 bohrs, and BEIN = 1.21
-
~
(49)
4
Stwalley, W. C.; de Llano, M. Z.Phys. D
1986, 2, 153.
4 4 4 5 5 5 5 6 6 7 7 7 8 8 8 8 8 8 9 9 9 10 10
10 10 10 10 11 11 11 12 12 12 13 13 13 13 14 19 79
-
4.82 4.43 4.82 4.55 4.47 linear 4.40 4.38 T-shaped rhombus 4.41 4.66 5.91 tetrahedron 4.40 4.82 4.45 linear 4.51 centered tetrahedron 4.24 sq pyram 4.40 4.82 4.60 trigonal bipyram dicapped 4.55 tetrahedron 4.47 octahedron 4.56 tricapped tetrahedron 4.54 fccl pentagonal bipyram 4.41 4.36 cubic 4.58 quad-capped tetrahedron square antiprism 4.41 fccl3 4.54 fccl7 4.54 fccl2 4.54 fcc136 4.55 fcc127 4.54 fcc123 4.54 fcc1368 4.56 fcc1357 4.55 fcc1238 4.55 fcc1236 4.55 fcc1235 4.55 4.54 fcc1234 fcc24578 4.56 fcc34568 4.55 fcc45678 4.55 fcc245678 4.57 fcc234568 4.55 fcc345678 4.56 fcc2345678 4.56 hcp structure 4.57 cubooctahedron 4.57 4.55 icosahedron 4.48 fcc12345678 4.58 capped 4.66 cubooctahedron Oh 4.82 4.80 4.82 bulk (fcc)'
60.0' 75.5O 81.8' 60.0' 4.33 bohrs 4.20 bohrs 4.44 bohrs 60.6' 55.8' 47.3'
4.30 bohrs 4.46 bohrs 4.41 bohrs
4.54 bohrs
1.06 0.89 0.52 0.39 1.21 1.23 1.34 1.47 1.14 0.60 1.56 1.18 1.26 1.25 1.35 1.67 1.57 1.71 1.80
DV-XLY 19'
1.84 1.88
CEM-N CEM-N
1.90 1.93 1.84 1.93
CEM-N CEM-N CEM-N CEM-N
1.94 1.95 1.96 1.97 1.99 2.00 2.01 2.01 2.04 2.04 2.05 2.05 2.06 2.07 2.08 2.09 2.11 2.11 2.12 2.14 2.22 2.22 2.19 2.28 2.17 2.38
CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N CEM-N
2.82 3.03
CEM-N
3.50
CEM-N
ECP-LC 7 ECP-CI 8
SCF
6
CEM-N
LSD
5
CEM-N CEM-N ECP-LC 7 ECP-CI 8 CEM-N
DV-XLY 19' LSD 5 CEM-N CEM-N CEM-N
DV-XLY 19' CEM-N CEM-N
LSD
5
CEM-N CEM-N CEM-N
LSD
5
'The u parameter for most geometries is defined in Figure 5 or Chart I of ref 1. * I f applicable, either the bond distance, b, or the bond angle, @, is defined in Figure 5. CBE/Nis the binding energy per atom. dThe experimental potential is predicted by definition, with the binding parameters obtained from ref 47. CThesecalculations did not optimize the atomic separations. fThe experimental potential is predicted by definition, with the binding parameters obtained from ref 44. The u parameter is the NND.
1564 The Journal of Physical Chemistry, Vol. 93, No. 4, 1989
Kress et al.
C u 5 D3h
-
\
5 80
J
?!.
t
510
J
E5
5.w
4: 3 4.60
140
ioc
1 '
1
L C
4
li
_"
5 bo
" ""
180
L 20
460
500
5 40
580
6 20
6 60
700
EQbATORIAL TO E 3 U A T O 3 A L ( a u )
Figure 9. Contour plot of BE for the variation of both the axial-toequatorial and the equatorial-to-equatorialdistances in Cu5. The arrows indicate the value of the axial-to-axial distance at the selected points.
eV, or slightly more stable than the linear triatomic system. The T-shaped structure is more stable with a B E I N = 1.34 eV at a = 4.38 bohrs and b = 4.44 bohrs. The more symmetrical rhombus structure yielded A = 4.41 bohrs, j3 = 60.6O, and B E I N = 1.47 eV. As in the Mg case, this structure is nearly a 60° rhombus which provides an additional nearest-neighbor ( N N ) interaction along its short diagonal in comparison to a square cluster. Allowing the four-atom planar structure to deform into a threedimensional tetrahedron increases the predicted B E I N to 1.56 eV, where the equilibrium N N D is found as 4.40 bohrs. This structure is the most stable of all the four-atom clusters. In comparison with other calculations, the ECP-LC method' yields B E I N = 1.14 eV for a rhombus, and the DV-Xcu approachI9yields B E I N = 1.I 8 eV for an unoptimized tetrahedron (with N N D = 4.82 bohrs, the bulk NND). The most extensive study to date is the local spin density (LSD) calculations of Delley et al.? where Cu clusters containing up to 79 atoms were considered, with many a t optimized geometries. For the linear structure, the LSD approach yields a = 4.40 bohrs, b = 4.20 bohrs, and B E I N = 1.23 eV. The tetrahedron calculated within the LSD approach is predicted to be more stable, as is the case within the C E M - N approach, with an equilibrium structure of N N D = 4.45 bohrs and B E I N = 1.26 eV. Although predicting the same trends, the LSD and CEM-N values for the B E I N are in disagreement by 0.3 eV. The most surprising feature is the radically different predictions of the relative stability of linear vs tetrahedral arrangements, 0.34 vs 0.03 eV for C E M - N and LSD, respectively. Moving on to the five-atom systems, for the linear structure we found a = 4.51 bohrs, b = 4.30 bohrs, and B E I N = 1.25 eV. As was the case for the five-atom Mg study, next we considered three-dimensional structures. For the centered tetrahedron, the N N D = 4.24 bohrs and B E I N = 1.35 eV, yielding only a small increase in stability. For the square pyramid, the equilibrium structure consists of a = 4.40 bohrs, b = 4.46 bohrs, and B E I N = 1.67 eV. The trigonal bipyramid is slightly more stable, as was the case for Mg, with the equilibrium values of a = 4.60 bohrs, b = 4.41 bohrs, and B E I N = 1.71 eV. There are no a b initio potentials of any accuracy to compare with mainly due to the sheer number of electrons contained in a five-atom Cu cluster. A HF-SCF calculation6 for this sysyem yields B E I N = 0.6 eV for either of the last two three-dimensional structures we have considered, which is clearly a severe underestimate. A nonoptimized DV-Xa c a l c ~ l a t i o nfor ~ ~the square pyramid, again at a = 4.82 bohrs, yields B E I N = 1.57 which is comparable to the predicted CEM-N result. As for the Mg cluster, the Cu5 Djhstructure is quite floppy as shown in Figure 9. In fact, the -7.9-eV contour allows a change of -3 bohrs in a. This is readily understood by noting that as the equatorial-to-equatorialdistance increases for a constant axial-to-equatorial distance, the axial-to-axial distance decreases. Consequently, the decreasing attraction between the equatorial atoms is partially compensated by the increasing attraction between the axial atoms.
Figure 10. Contour plot of BE for the variation of both the axial-toequatorial and the equatorial-to-equatorialdistances in Cu,. The arrows indicate the value of the axial-to-axial distance at the selected points.
There are no high-quality geometrically optimized calculations for CuN, N = 6-12, and thus we have spent considerable effort on these structures to provide such information. On the basis of our previous experience, we have considered only three-dimensional structures with large average coordinations of their atoms. One point of notation is worth emphasizing. The octahedron is the basic building block of many larger structures formed by placing atoms at the corners of the cube circumscribing the octahedron. These structures are referred to as fccpqrs ..., wherep, q, etc., are the locations of the corner atoms. The results are listed in Table 111. For the six-atom system the equilibrium values for the dicapped tetrahedron were N N D = 4.55 bohrs and B E I N = 1.80 eV, while those for the octahedron were N N D = 4.47 bohrs and B E I N = 1.84 eV. A similar decrease in the N N D and increase in the B E I N between the C , and oh structures was observed for the Mg clusters. Although CEM-N has predicted generally that the geometry with the largest average coordination as the most stable for both Cu and Mg clusters, these two geometries provide the same average coordination of four. The slight difference in energy between these two six-atom structures is attributed to differences in second N N interactions: the ratio of the second N N distance to the N N D is'larger for the C, geometry than for the Oh geometry. This second N N effect is small since the difference in B E I N is only 0.04 eV. The D5h seven-atom cluster provides an interesting contrast to the D3h five-atom cluster. As shown in Figure 10 the equatorial-to-equatorial distance cannot increase very much, keeping the axial-to-equatorial distance constant, before the axial-to-axial separation becomes so small as to lead to a repulsive interaction between the axial atoms. The contrast with the trigonal bipyramid in Figure 9 is quite startling and arises from the decreased axial-to-axial separation near the equilibrium for the pentagonal bipyramid as indicated in the two figures. The eight-atom structures are interesting since B E I N varies only slightly while the N N D changes considerably among the different clusters. This feature does not recur for any of the larger structures that we treated, but these are all rather closely packed, being mainly additions to the octahedron. The 13-atom cubooctahedron was the subject of a high-quality LSD calculation? and thus we single it out for special discussion. The potential was calculated as a function of the NND, a so-called spherical breathing mode, and is shown in Figure 11. The minimum is at N N D = 4.57 bohrs with B E I N = 2.22 eV. Note how flat the potential is between the bulk distance at 4.82 bohrs and the minimum at 4.57 bohrs and note the excellent agreement between the C E M - N and LSD methods. The major advantage of the C E M - N method is that it is much cheaper in terms of computational effort, since the LSD method is an all-electron S C F type of calculation. A typical calculation of the energy in the 13-atom system without using any symmetry (to reduce the time required to evaluate AG) requires about 32 min on a Ridge 3200 (about 3 times faster than a V A X 11/780), while the all-electron
The Journal of Physical Chemistry, Vol. 93, No. 4, 1989 1565
Corrected Effective Medium Method
CEM-N CU POTENTIAL h
/
CEM-N ( 4 . 5 7 . 2 . 2 2 )
/
2.3 m
BULK LRTTICE
4.2
4.8
4.4
I1
4.8
MRREST NEIGHBOR DISTANCE (bohr)
Figure 11. Variation of BEIN with NND for the Cu13cubooctahedron with Ohsymmetry. The value of BEIN for a NND distance equal to the bulk value is shown by the vertical arrow. The asterisk indicates the geometry-optimized point from ref 5 . The HF-SCFvalue of ref 6 is only noted, since the BEIN is off scale, but it does yield a NND in agreement with the other calculations.
-
Cu,, C I us te rs 4.0
> v a,
BULK
201 0
z 3.0 h
t-
N
(22
LL E
z
w
s
135
\
=
430
19’ ;I3
2.0
(22
z,
2
1.0
03 v
0.0 0.0
o
CEM-N
(bulk d i s t . )
x
LSD
(optim. dist.)
0.4
0.2
X
=
1
-
(1
0.6
- N-”3
0.8
1.0
)3
Figure 12. CEM-N binding energy per atom for N-atom Cu clusters constructed as fixed “spherical” fragments of the bulk fcc lattice (see text) for N 2 13 and geometrically optimized for N = 2-13. All calculations used the covalent embedding function. The values are plotted as a function of the fractional surface volume, X. Also provided are the values from local spin density (LSD) calculations by Delley et aL5
LSD results must utilize significant time on a class VI supercomputer. With a nonoptimized implementation of the better integration scheme of ref 43b this time has now been reduced to about 4 min. Furthermore, the C E M calculations increase in computational time only slightly faster than linearly in the number of atoms as long as the evaluation of AG dominates the computational requirements. At this point we would like to point out a couple of trends in Table 111 for the Cu clusters. All of the N N D s for the clusters
considered fall between a lower limit, that of the diatomic bond distance, and an upper limit, that of the bulk metal NND. The general trend is for an increasing N N D with cluster size. Also note that the B E I N increases as function of N. Thus, the properties of the Cu clusters vary much more smoothly as a function of N than these properties did for the Mg clusters. To examine the trends for the binding potentials for larger Cu clusters, we also have considered “spherical” fragments of the bulk fcc lattice, where we have constrained the cluster N N D to be equal to the bulk NND. This is directly analogous to the large cluster calculations for Mg. The values of B E I N that we have calculated with the CEM-N method are plotted in Figure 11 as the same function of N we considered for the large Mg clusters (i.e., X , the fractional surface volume 49). The CEM-N values in Figure 11 for N I13 correspond to the optimized values for the most stable structures in Table 111, with the other C E M - N values corresponding to the fixed “spherical” fragments. As in the Mg case, a (linear) extrapolation to X = 0 ( N m) of the CEM-N cluster data (-4.1 eV) overestimates the actual bulk cohesive energy by about 0.6 eV. This again suggests that the X variable may not be the one of choice for this extrapolation. In Figure 11 we have also provided the optimized LSD result of Delley et al. for N = 13, as well as their N = 79 result for a nonoptimized cluster with a N N D fixed at 4.80 bohrs, which yields B E I N = 3.03 eV. The agreement between the N = 79 values in Figure 11 is excellent.
-
IV. Summary and Conclusions We have presented the corrected effective medium theory formalism for the calculation of the energy per atom of an infinite 3-D periodic system. The determination of the embedding functions for homogeneous systems was facilitated by using both diatomic binding potentials and bulk cohesive energies. The two different systems determined low- and high-density parts of the covalent embedding function, respectively, with the two parts smoothly connected via interpolation. This provided strong evidence for the universality of the embedding energy function as a function of the size of the system, at least for the homogeneous atomic systems considered here. Using these covalent embedding functions, we investigated the structures and energies of a number of Mg and Cu clusters containing from 2 to 201 atoms. Potential energy surfaces for isomerization of various small clusters were determined, and some clusters were shown to be especially “floppy”. The equilibrium geometries of moderate-size clusters were shown to be significantly contracted with respect to the bulk distance. This feature could be of critical importance in the description of reactions between gases and small metal clusters, in which the analogous gassurface reaction is found to be structure sensitive. Acknowledgment. This work was supported by the Director for Energy Research, Office of Basic Energy Sciences. Ames Laboratory is operated for the US. Department of Energy by Iowa State University under Contract No. W-7405-ENG-82. Registry No. Cu, 7440-50-8; Mg, 7439-95-4.
