© Copyright 1996 by the American Chemical Society
VOLUME 100, NUMBER 26, JUNE 27, 1996
LETTERS Local, Gradient-Corrected, and Hybrid Density Functional Calculations on Pdn Clusters for n ) 1-6 Gabriele Valerio and Herve´ Toulhoat* Institut Francais du Petrole, 1 et 4 AVenue de Bois Pre´ au, BP 311 92506 Rueil-Malmaison, France ReceiVed: February 8, 1996; In Final Form: April 30, 1996X
Full optimizations of Pdn cluster geometries have been performed using density functional methods and allowing the occurrence of low-symmetry structures. To account for at least partially for relativistic effects, use was made of Hay and Wadt RECP. The following ground states were found respectively for n ) 1, 2, 3, 4, and 6: 1S(4d10), 3Σ+u , 3B2 {C2V}, 3B2 {C2V}, 3B2g {D4h}. Various exchange-correlation functionals have been tested on Pd2 and Pd3, including the recent B3LYP one. The introduction of some exact exchange reduces significantly the predicted binding energy, which thus matches the experimental values in the case of Pd2. Moreover, it allows us to overcome the serious SCF convergence problems arising from the presence in these clusters of a large number of low-lying, nearly degenerate states.
1. Introduction Transition-metal clusters are the topics of several experimental and theoretical studies.1-4 Small clusters are suitable models to simulate surface properties and chemisorption. This is particularly true when the reactivity depends on local features of the geometrical structure and when the real periodicity of the surface does not influence the local electronic properties. On the other hand, cluster studies are interesting to understand the properties of small metallic aggregates, formed on catalytic supports or in the gas phase. From a more fundamental point of view, the theoretical electronic structure investigations on small clusters (even dimers) give important information about the nature of metal-metal bond and cluster-adsorbate interactions. Experimentally, supported catalysts preparation techniques on one side and condensation of atomic metal vapor produced by pulsed laser on the other side yield very small particles for which many properties can now be measured and calculated. Properties such as ionization potential, electron affinity, reactivity,5-7 and magnetic behavior8,9 are studied X
Abstract published in AdVance ACS Abstracts, June 15, 1996.
S0022-3654(96)00356-5 CCC: $12.00
experimentally on these clusters. Cluster sizes having a particular stability (magic numbers) can also be observed using mass spectroscopy.2 On the contrary, the geometrical structure is difficult to determine experimentally, especially in the gas phase; in this case the theoretical investigation becomes the better tool to know the equilibrium geometry of metallic clusters. Clusters often present properties that are strongly sizedependent and that bear no resemblance to the crystalline state. To investigate these properties and their size dependance, by means of electronic structure calculation methods, the most stable structure must be determined. We need to relax completely all the geometrical degrees of freedom, as well as to determine the electronic ground state. This is accomplished looking for the true minimum of the cluster potential energy surface; this surface is a multidimensional one and usually includes several minima. A very rigorous investigation can be carried out connecting the molecular dynamic methods to the ab initio electronic structure calculations10,11 and making geometrical optimization by means of the simulated annealing technique. This kind of investigation is moreover computationally very expensive. It is preferable to use first models with © 1996 American Chemical Society
10828 J. Phys. Chem., Vol. 100, No. 26, 1996 a lower level of computational intensity. More accurate approaches must then be used in order to detect small structural distortions (e.g., due to Jahn-Teller effects).12 Stave and DePristo13 have applied CEM (corrected effective medium) methods14 to palladium and nickel clusters up to 23 atoms. Their general structural findings are high symmetry of stable clusters (e.g., Td, Oh, Ih, and Td for 4, 6, 13, and 17 atoms, respectively) and shorter bond distances compared to crystal. They found that structural properties of transition-metal clusters can be very different from that of the bulk and are often characterized by 5-fold symmetry, which is forbidden in crystals, owing to the translational symmetry. These results are a good starting point for more accurate calculations. The electronic ground-state determination of transition-metal clusters is particularly difficult, because of the abundance of low-lying electronic states and isomers due to incomplete d-shells. Electronic correlation and relativistic effects play a fundamental role. For this reason there are few high-level ab initio studies on transition-metal clusters that include complete geometry optimization. In particular, in the case of palladium, the only ab initio investigation on Pd3 and Pd4, where the possibility of lower symmetry is considered, are those due to Balasubramanian,15,16 even if dimer was studied at postHartree-Fock level and density functional (DFT) level. Some preliminary results, which show marked differences with respect to those reported in ref 15, are also obtained by Blomberg et al.17 Fahmi and van Santen18 have made an interesting study of ethylene adsorption on palladium clusters of different sizes, relaxing atomic positions but fixing the higher symmetry for bare clusters (D3h for Pd3, for example). Therefore, we have decided to start a characterization of Pdn (n ) 1, 2, 3, 4, and 6) with the aim of testing the calculation methods and comparing the results of our DFT study with the HF/CI study15 on Pd3. We have also investigated the occurrence of low-symmetry structures for n ) 4, 6. 2. Computational Details The calculations reported here were carried out mostly using the density functional options of Gaussian94.19 To incorporate the scalar relativistic effects, we have used the LanL2DZ19 basis set (Hay and Wadt relativistic effective core potentials (RECP), small core20 plus DZ). Geometry optimizations were performed at local (L) and nonlocal (NL) levels, for comparison. Gaussian94 offers a wide variety of DFT models, LSDA and GGA, as well as hybrid methods,21,22 which are characterized by a mixing of exact exchange (“Hartree-Fock”) and exchangecorrelation from DFT. The straightforward substitution of DFT approximate exchange by exact HF exchange has proven to be inadequate for molecules.23 To consider the exact exchange information, Becke has proposed a hybrid method21 with a new partitioning of exchange and correlation energies between HF and DFT terms. This partitioning is based on an approximation of the rigorous “adiabatic connection” formula, for the exchangecorrelation energy in the Kohn-Sham theory. The following functional was defined and parametrized by Becke:22 LSDA EXC ) EXC + a0(EXexact - EXLSDA) + aX∆EXB88 + aC∆ECPW91 (1)
where a0 ) 0.20, aX ) 0.72, aC ) 0.81, EXexact is the exact exchange energy, ∆EXB88 is Becke’s gradient correction for is the gradient correction for correlaexchange,24 and ∆EPW91 C tion of Perdew and Wang.25 The Becke three-parameter functional (B3LYP26 in Gaussian94) incorporates 20% of exact
Letters TABLE 1: Results for the Palladium Dimer Ground-State Triplet and Lower Excited-State Singleta triplet
singlet
method
Re
De
Re
De
CASSCF/MRSDCI15,28 SCF/MCPF17 DFT, L/NL29 DFT, NL/NL18 Dmol, L/NL G94, DFT, L/NL G94, DFT, NL/NL G94, hybrid, NL/NL experimental34,36
2.48 2.55 2.46 2.56 2.48 2.44 2.54 2.53
0.86 0.47-0.64 1.35 1.27 0.05 1.31 1.36 0.96 0.73-1.13
2.87 2.88 2.64 2.74 2.69 2.61 2.74 2.76
0.69 0.16 0.95 1.09 0.37 0.96 1.00 0.62
a R is the equilibrium bond distance (Å); D (eV) is calculated with e e respect to Pd 1S(d10) atoms. L ) local, NL ) nonlocal, approximation. L/NL means the geometry was optimized at local level, but NL singlepoint calculations were made to obtain De.
exchange in the GGA functional. We have applied the Becke method to small palladium clusters. The electronic structure of palladium clusters is characterized by a large number of low-lying, nearly degenerate states. This involves several problems when quantum calculations are made on these systems. First, the true electronic ground state is very challenging to find. Actually it is easy to fall in local minima, which correspond to geometry very different from that of the cluster in the ground state. Second, when the HOMO-LUMO gap is small, SCF convergence problems arise. This is particularly true in the case of DFT methods, which tend to underestimate this gap. The hybrid methods, that yield more reasonable HOMO-LUMO energy differences, are an operating way in the investigation of transition-metal cluster properties, so we have decided to consider Pd2 as a test for the performances of Becke functional.22 The results show that the SCF convergence problems are removed. The triplet-singlet relative stability is well reproduced, with an energy separation similar to the GGA (BLYP) value. In the same way, the calculated energy difference between the 1S(4d10) ground state and the 3D(5s14d9) excited state of the palladium single atom is 0.87 eV, either using the B3LYP functional or the GGA (BLYP) one. The experimental value is 0.82 eV.27 The interatomic distance for Pd2 ground state calculated using this functional is longer compared to our LSDA results and to accurate results recently reported in the literature12,28,29 but compare well with MCPF Blomberg results17 and with GGA geometries (see Table 1 and refs 18 and 30). Some results obtained with the Dmol program,32 which employs localized numerical orbitals, are also reported for comparison and in order to display the importance of relativistic effects in palladium compounds. 3. Results and Discussions 3.1. Pd2. The dimer represents the smallest system which allows us to investigate the metallic bond. Owing to its small size, the dimer allows one to exploit the most accurate theoretical methods, as well as to perform easily and accurately geometry optimizations. Pd2 is then a good system to test the methods we want to use on bigger clusters. Recent Pd2 electronic structure calculations12,17,18,28-30 have all established a triplet ground state (3Σ+ u ) and a singlet lowest ). The results we obtain with Dmol32 give an excited state (1Σ+ g inverse stability, with a lower energy for the singlet, even though the respective geometries are correct. Binding energies (Table 1) are weak, compared to experimental data. The Mulliken
Letters
J. Phys. Chem., Vol. 100, No. 26, 1996 10829
TABLE 2: Mulliken Population Analysis for the Triplet and Singlet Electronic States of Pd2
TABLE 4: Energy Separations (kcal/mol) of Triplet Electronic States of Pd3a
5s/5p/4d method
triplet
singlet
CASSCF/MRSDCI28 DFT29 DFT18 Dmol G94, DFT G94, hybrid
0.65/0.07/9.28 0.58/0.04/9.38 0.55/0.02/9.43 0.53/0.05/9.43 0.56/0.07/9.36 0.58/0.05/9.36
0.10/0.05/9.86 0.17/0.02/9.80 0.11/0.04/9.85 0.26/0.07/9.69 0.13/0.07/9.81
TABLE 3: Results for Pd2, Obtained Using the Standard or Hybrid DFT Options of Gaussian94a A B C
triplet singlet T-S (eV) triplet singlet T-S (eV) Triplet singlet T-S (eV)
Pd-Pd
Pd2
Pd
∆E
2.44 2.61
-253.2758 -253.2629 -0.351 -253.2774 -253.2642 -0.360 -253.4488 -253.4362 -0.343
-126.5817 -126.6138 0.874 -126.5817 -126.6138 0.874 -126.6749 -126.7068 0.866
-1.31 -0.96
2.54 2.74 2.53 2.76
BLYP B3LYP ref 15
populations analysis (Table 2) shows a significant population of 5s orbitals in the triplet, differently from the singlet case. The analysis done by Balasubramanian28 proves that most of the electronic states arise from d9s1-d10 mixture of atomic configurations, with d9s1 predominance; an exception to this is the 1Σg+ state which was found to be predominantly arising from d10 configuration. Therefore, the wrong results obtained with Dmol seem to derive from relativistic effects not being accounted for. We can expect in reality a 5s contraction resulting from the massvelocity correction would favor the binding between palladium atoms through the superposition of 5s orbitals and induce the stabilization of the triplet with respect to the singlet state. The contraction of the outer valence 5s orbital of the palladium atom hence induces enhanced stability and shorter Pd-Pd bond, like in the case of gold compounds.33 The same findings were reported by Dixon and co-workers.30 The calculations made using relativistic ECP within the DFT options of Gaussian94 give very good results. In particular those obtained with the hybrid functional are very encouraging. The B3LYP interatomic distances (Tables 1 and 3) are similar to the BLYP and MCPF17 results. In addition we have less convergence problems, a better HOMO-LUMO separation, and therefore large time saving for SCF convergence. Moreover we have to notice that the Pd-Pd interaction energy calculated with the B3LYP functional is lower than the GGA one (Table 3) and close to the experimental values (Table 1). 3.2. Pd3. The triangular structure of Pd3 has been fully optimized in the C2V symmetry. Previous investigations have pointed out the triangular geometry is more stable for Pd3 relative to the linear geometry15,18 and the symmetry could be lower than D3h.15 We have performed geometry optimization using the BLYP and B3LYP functionals. Even for the Pd3 case, the hybrid functional makes SCF convergence and ground-state determination easier with respect to standard DFT calculations. In either cases geometry results are similar, while the B3LYP atomization energy (2.50 eV) is quite smaller than the BLYP value (3.66 eV). The relative energies between the different
0.5 0 0
2.8 2.0 0.5
3A
1
2
3B 1
0.8 3.0 1.5
0 3.8 4.5
TABLE 5: Geometries and Total Energies for the Electronic States of Pd3, Obtained Using the B3LYP Functionala
a
a A: LSDA(VWN) geometries and GGA (BLYP) energies. B: BLYP geometries and energies. C: B3LYP geometries and energies. ∆E (eV) is calculated with respect to Pd 1S(d10) atoms. Total energies are in hartree, distances in angstroms.
3A
a Comparison between DFT (BLYP), hybrid (B3LYP), and postHartree-Fock results.
-1.36 -1.00 -0.96 -0.62
3B 2
state
R
θ
E
3B 2 3A 1 3A 2 3 B1 1 A1
2.545 2.612 2.575 2.567 2.516
65.6 57.8 60.0 60.0 60.0
-380.2144 -380.2111 -380.2096 -380.2083 -380.2048
Total energies are in hartrees, distances (R) are in Å.
triplet states we have found have proved a nearly inverse stability (Table 4). The energetic order given by the hybrid functionals for the triplet states of Pd3 follows closely the one obtained by Balasubramanian,15 whose work is the only extensive study on the Pd3 cluster. He has investigated, with the MRSDCI method and LaJohn RECP,35 several low-lying electronic states of Pd3 in the C2V symmetry. The main difference between his results and ours is that he has found a singlet (1A2) ground state, with nonequilateral geometry (Re ) 2.67 Å; θe ) 55°); the lower excited state he has found is the triplet 3B2, 6.5 kcal/mol higher than the GS and with a geometry (Re ) 2.52 Å; θe ) 67.5°) similar to that we have obtained for this state (see Table 5). Blomberg17 post-Hartree-Fock calculations resulted in a 3B2 GS, in contrast with ref 15 and in agreement with our results using the hybrid functional. The lower singlet state we determine (Table 5) has D3h (1A1) symmetry and is 5.6 kcal/mol higher than the 3B2 GS. The interatomic distances we have determined within the local approximation are always shorter than those of ref 15. The role of 5s orbitals in the Pd-Pd bond is significant, for all the electronic states we have studied. The 5s population is higher the shorter the Pd-Pd distance. Therefore, the 5s population is always higher for the apex atom (Pd1), except for 3A1 which has the only geometry with θ < 60° and so a shorter Pd2-Pd3 bond. 3.3. Pd4 and Pd6. All the calculations for these two clusters were performed with the B3LYP functional. The full geometry optimization starting from the tetrahedral Pd4 gave a C2V JahnTeller distorted structure (3B2 ground state), in accordance with ZINDO calculations by Estiu` and Zerner.12 The six bond lengths are respectively 2.60, 2.78, and 2.64 Å (×4). The lower excited state we have found is a 1A1 state with S4 point group (Pd-Pd ) 2.85 Å (×2) and 2.58 Å (×4)), 0.70 eV higher in energy than the GS. Even for Pd6 the cluster has not the full Oh symmetry. In this case we have determined a triplet GS (3B2g) with D4h symmetry (one axis of the octahedron is longer than the other two) and the following bond lengths: 2.64 Å (×4) and 2.75 Å (×8). The first singlet state (1Ag) has the same structural symmetry (Pd-Pd ) 2.58 and 2.81 Å) and is 0.95 eV higher than the GS. Table 6 gives the calculated atomization and binding energies, with other properties for all the clusters studied. The energy differences are relative to atomic palladium in its ground state 1S(4d10). Even if the atomization energies increase with increasing atom number, the binding energies decrease and fall between Pd2 and bulk experimental values.36
10830 J. Phys. Chem., Vol. 100, No. 26, 1996 TABLE 6: Atomization Energies (AE), Pd-Pd Bond Energies (BE), and Averaged Mulliken Populations, Obtained Using the B3LYP Hybrid Functionala
Pd2 Pd3 Pd4 Pd6 bulk
m
Pd-Pd
AE
BE
1 3 6 12
2.53 2.55-2.76 2.60-2.78 2.64-2.75 2.75
22.1 59.0 114.9 195.2 90.0
22.1 19.7 19.2 16.2 15.0
Mulliken populations 5s/5p/4d 0.58/0.05/9.36 0.62/0.08/9.22 0.59/0.11/9.30 0.47/0.17/9.36
a AE ) E[Pd ] - nE[Pd]; BE ) AE/m. m is the number of Pd-Pd n bonds in the cluster. Bulk data taken from the Handbook of Chem. Phys., 1995. The bond lengths are in angstroms, energies in kcal/mol.
This is in agreement with the statement of the bond order conservation principle37 and with the increased connectivity of each Pd atom. Similar results are found by Fahmi and van Santen.18 On the contrary, Balasubramanian has found an impressive enhancement of stability from Pd2 to Pd3 and so a very high Pd-Pd bond energy for the latter, even higher than for Pd4.16 We note also that mean bond lengths are shorter in these clusters than in the bulk (2.75 Å). The averaged Mulliken populations (for n ) 3, 4, 6) display an increasing p and d population and a decreasing s, according to the lengthening of Pd-Pd distances and the growth of the coordination number. 4. Conclusions Research on cluster catalytic activity of palladium and transition metals in general needs the knowledge of bare clusters ground-state properties. High-level first principles investigations are scarce, because of calculation costs, but also due to the SCF convergence and ground-state determination difficulties. Becke’s hybrid methods seem promising tools to overcome these problems, within the density functional formalism. We have applied them to Pdn clusters (n ) 1, 2, 3, 4, 6), obtaining quite encouraging results. The binding energies obtained using the Becke three parameters hybrid functional (B3LYP) are about 30% smaller than those calculated with standard gradient corrected functionals. They are still greater than the bulk cohesive energy, but quickely approaching it. For n ) 3, 4, 6 structural distortions from perfect symmetries were found (respectively C2V, C2V, and D4h instead of D3h, Td, and Oh), presumably due to Jahn-Teller effects. Moreover, we have verified that not only electronic correlation but also relativistic effects have to be considered in every first principles study on palladium clusters, because of the close energies between (n + 1)s and nd atomic orbitals. Acknowledgment. The authors would like to thank Dr. E. Wimmer and Dr. C. Guerra, as well as Professor G. Pacchioni for helpful discussions. Professor R. A. van Santen and Dr. A. Fahmi are kindly acknowledged for a preprint of ref 18. G.V. thanks the Institut Franc¸ ais du Pe´trole for a postdoctoral fellowship.
Letters References and Notes (1) Pacchioni, G.; Bagus, P. S.; Parmigiani, F. Cluster Models for Surface and Bulk Phenomena, NATO ASI series B, Physics; Plenum: New York, 1992; Vol. 283. (2) Koutecky, J.; Fantucci, P. Chem. ReV. 1986, 86, 539. (3) Morse, M. Chem. ReV. 1986, 86, 1049. (4) Schmid, G. Chem. ReV. 1992, 92, 1709. (5) Ervin, K. M.; Ho, J.; Lineberger, W. C. J. Chem. Phys. 1988, 89, 4514. (6) Fayet, P.; Kaldor, A.; Cox, D. M. J. Chem. Phys. 1990, 92, 254. (7) Ren, X.; Hintz, P. A.; Ervin, K. M. J. Chem. Phys. 1993, 99, 3575. (8) Milani, P.; de Heer, W.; Chatelaˆin, A. In p 67 of ref 1. (9) Douglass, D. C.; Bucher, J. P.; Bloomfield, L. A. Phys. ReV. B 1992, 45, 6341. (10) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (11) Car, R.; Parrinello, M. Phys. ReV. Lett. 1988, 60, 271. (12) Estiu´, G. L.; Zerner, M. C. J. Phys. Chem. 1994, 98, 4793. (13) Stave, M. S.; DePristo, A. E. J. Chem. Phys. 1992, 97, 3386. (14) Kress, J. D.; Stave, M. S.; DePristo, A. E. J. Phys. Chem. 1989, 93, 1556. (15) Balasubramanian, K. J. Chem. Phys. 1989, 91, 307. (16) Dai, D.; Balasubramanian, K. J. Chem. Phys. 1995, 103, 648. (17) Blomberg, M. R. A.; Siegbahn, P. E. M.; Svensson, M. J. Phys. Chem. 1992, 96, 5783. (18) Fahmi, A.; van Santen, R. A., preprint. (19) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; HeadGordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, Revision B.3; Gaussian, Inc.: Pittsburgh, PA, 1995. (20) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299. (21) Becke, A. D. J. Chem. Phys. 1993, 98, 1372. (22) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (23) Ziegler, T. Chem. ReV. 1991, 91, 651. (24) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (25) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. ReV. B 1992, 46, 6671. (26) In the B3LYP option, the gradient corrections for correlation of Lee-Yang-Parr are used, instead of that of Perdew-Wang. Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (27) Moore, C. E. Table of Atomic Energy LeVels; U.S. National Bureau of Standards, Washington, DC, 1971. (28) Balasubramanian, K. J. Chem. Phys. 1988, 89, 6310. (29) Goursot, A.; Papai, I.; Salahub, D. R. J. Am. Chem. Soc. 1992, 114, 7452. (30) Nakao, T.; Dixon, D. A.; Chen, H. J. Phys. Chem. 1993, 97, 12665. (31) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (32) Delley, B. J. Chem. Phys. 1990, 92, 508. Dmol is available from Biosym Technologies, San Diego. The 2.36 version was used, with Double Numerical basis set and FINE parameter for the mesh. (33) Balasubramanian, K. J. Phys. Chem. 1989, 93, 6585. Hay, P. J.; Wadt, W. R.; Kahn, L. R.; Bobrowicz, F. W. J. Chem. Phys. 1978, 69, 984. (34) Lin, S. S.; Strauss, B.; Kant, A. J. Chem. Phys. 1969, 51, 2282. (35) LaJohn, L.; Christiansen, P. A.; Ross, R. B.; Atashroo, T.; Ermler, W. C. J. Chem. Phys. 1987, 2812. (36) Handbook of Chemistry and Physics, 76th ed.; Lide, D. R., Ed.; CRC Press: Boca Raton, FL, 1995. (37) Shustorovich, E. Surf. Sci. Rep. 1986, 6, 1.
JP960356Q
