27 June 2002
Chemical Physics Letters 359 (2002) 428–433 www.elsevier.com/locate/cplett
The importance of global minimization and adequate theoretical tools for cluster optimization: the Ni6 cluster case Fernando Ruette b
a,*
, Carlos Gonz alez
b
a Laboratorio de Quımica Computacional, Centro de Quımica IVIC, Apartado 21827, Caracas, Venezuela Physical and Chemical Properties Division, NIST, Room A11, Building 211, Gaithersburg, MD 20899, USA
Received 20 December 2001; in final form 5 April 2002
Abstract A Monte Carlo simulated annealing (MCSA) algorithm using parallel processing has been employed in quantum parametric method in order to full optimize the Ni6 cluster. Results obtained from this procedure show that the same trend is reproduced by standard DFT optimization calculations. However, results do not correspond with the octahedral geometry reported in several theoretical DFT calculations and for classical global minimization techniques. Different structures close in energy are found here, as the most stable ones: pentagonal pyramid (PPY) and capped trigonal bipyramid (CTB). Ó 2002 Elsevier Science B.V. All rights reserved.
The interest for small metallic clusters in technological issues may be assessed by the number of recent publications. For example, dozens of papers related with nickel metallic clusters [1,2] have been published in the last three years. Applications of transition metal cluster chemistry have an enormous importance in different technological disciplines. It comprises synthesis of catalysts, in which size-controlled transition metal clusters are supported on porous materials. Specific electric and ferromagnetic properties of small metallic aggregates become of great importance in recording industry and high technology nanoelectronic de-
*
Corresponding author. Fax: +58-2-5041350. E-mail address:
[email protected] (F. Ruette).
vices such as: quantum dot lasers, single electron transistors, and diminutive magnets. In many cases, accurate electronic and structural information of metallic clusters is impossible to obtain from experimental techniques because of the cluster size, and the fact that in some cases those are embedded into porous materials (catalysts). Therefore theoretical modeling is one of the most adequate tools to tackle this type of problem. Nevertheless, there is a compromise between method accuracy and model feasibility, because the selection of the proper theoretical method depends on the size of representative model. One feature that is closely related with most of the electronic properties of materials is the geometrical structure. In most of the cases for systems that contain transition metals, such as clusters,
0009-2614/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 ( 0 2 ) 0 0 7 3 3 - 9
F. Ruette, C. Gonzalez / Chemical Physics Letters 359 (2002) 428–433
full-geometry optimization for well-defined ab initio methods is prohibitive, because of manylocal-minima existence. The most standard methodology for searching the potential energy surface is based on molecular mechanics using empirical potentials. Even with this very simplified approach, in some cases, it can be very difficult task to find the global minimum. Reported geometries for metallic cluster can be very controversial using different type of methods. For example, several authors have assumed that Ni6 is octahedral [3–5] based on qualitative bulk information and from chemical probe techniques [6]. On the other hand theoretical approaches predict Oh [7–12], C5v [13], and D4h (distorted octahedral) [14] symmetries. In this letter we discussed a methodology to calculate global minimum and its application to the Ni6 case. A comparison of the most stable geometries of the Ni6 cluster was carried with different theoretical approaches: parametric quantum and DFT methods. Here, we will explore a relatively simple methodology for searching different geometric structures that for practical purposes at a higher level of theory it is very arduous to obtain the global minimum. Quantum parametric methods are based on the modeling of parametric Hamiltonians (Hpa ) for atoms and molecules. The Hpa is optimized with respect to one considered as exact (Hexa ), using simulation techniques. In previous works [15–17], it was shown that optimization of Hpa for molecular systems can be performed in terms of binding energies of X–Y bonds (BEXY ) from experimental data or results of very well-defined methods. X and Y stand for molecular fragments, as well as for atoms that form the X–Y bond. XY XY min Hexa Hpa ! 2 1=2 X XY XY ¼ min ; ð1Þ BEexaI BEpaI I
where XY XY X X X H W W H W BEXY ¼ WXY I I Y Y Y W H W :
ð2Þ
429
This formalism has been implemented in a modular code for studying surface-adsorbate interactions called CATIVIC [18], based on MINDO/SR [19]. Here, several parametric functionals [20], new computational tools, and different parametrization methods were included. Ni–Ni parameterization was performed with respect to new experimental values of the diatomic molecule [21]. The molecular parameters were obtained after the application of the above scheme, Eqs. (1) and (2). The calculated values were b ¼ 0:996204, for electronic terms and a ¼ 0:990750 and CNi–Ni ¼ 0:649851, for the core–core potential, as expressed in [16]. Results show that the most stable configuration corresponds to the r2g p4u d3g d3u p4g r2u r2g molecular orbital occupation with a multiplicity equal to 3. This state is the most stable experimentally reported, with a Ni–Ni equilibrium distance of and a binding energy of )47.0 kcal/mol 2.15 A [21]. The binding energy between metallic atoms is calculated according to the diatomic binding energy (DBE) obtained from partition total binding energy, as described elsewhere [22]. Monte Carlo simulated annealing (MCSA) [23] has been applied to a variety of problems related with quantum chemistry [24] and it is quite popular in global optimization and molecular dynamics studies. It is a technique based on stochastic search and the Metropolis algorithm [25]. In this work a new version of a parallel seriallike algorithm based on the Roussel–Ragot– Dreyfus method [26] (RRD) is used. This method is a hybrid type one, because high and low temperature modes can be merged onto one single algorithm. The advantage of this new algorithm lays on the improvement of the performance because the amount of communication between ‘master’ and ‘slaves’ that appear in the RRD method is significantly diminished. Each processor contains the information of the master processor about a starting configuration of coordinates and evaluates a new cost function with the corresponding Metropolis test. At high temperature regime, a new configuration is chosen at random from all acceptable moves obtained from different processors. On the other hand, at low temperature regime, the first accepted move is taken the re-
430
F. Ruette, C. Gonzalez / Chemical Physics Letters 359 (2002) 428–433
sulting geometry distributed to the rest of the processors [27]. MCSA was implemented on a massively parallel machine CRAY-T3D with 512 processors. All calculations were started at an initial temperature of 2000 K and quenched to 0.5 K using a fixed maximum step of quenching rate of 0.99 A 0.1 A was used in all simulations. The ‘initial guess’ was determined at random, and for a particular cluster size, a total of 100 runs were carried out. On each run, a total of 2000 equilibration steps were taken for each temperature. This number of steps was found to comply with the ‘equilibrium’ criterion previously discussed. This procedure was carried out for calculations with parallelized MCSA code using the parametric method CATIVIC.
CATIVIC calculations were performed with different multiplicities (5, 7, and 9). Results indicate that the optimal multiplicity was 7 and the global minimum corresponds to a pentagonal pyramid (PPY) geometry, as shown in Fig. 1a. This multiplicity is in agreement with other theoretical calculations [8,29]. Calculations for three different geometries (Figs. 1a and c were carried out using DFT approach from the Gaussian code [28] with corrected gradient and different basis sets: B3LYP/ LANL2DZ and B3LYP/LANL2MB. Values of energy differences for these clusters are presented in Table 1. These differences change with the basis set quality and the stability order of isomers is similar only for parametric and B3LYP/ LANL2MB calculations. Results suggest that the
Fig. 1. Geometries of Ni6 clusters: (a) pentagonal pyramid (PPY); (b) octahedral (OCT); (c) capped trigonal bipyramid (CTB).
F. Ruette, C. Gonzalez / Chemical Physics Letters 359 (2002) 428–433
431
Table 1 Total energy differences for Ni6 with OCT, CTP, and PPY structures Method structure
CATIVIC
DFT-ECP B3LYP-LANL2MB
DFT-ECP B3LYP-LANL2DZ
OCT CTP PPY
0.00 )48.8 )56.8
0.00 )6.69 )18.63
0.00 )11.29 )10.77
Values are in kcal/mol.
most stable geometry is not the OCT, as has been normally accepted [3–12]. It is expected that some classical potentials are bias to OCT structure because they, in general, are adjusted to mimic bulk properties. Our results of B3LYP/LANL2DZ and CATIVIC are in agreement with ab initio calculations of Nygren et al. [13] using one-electron effective core potential (ECP) and core polarization potential (CPP) with valence correlation through average coupled pair functional (ACPF). They conclude that PPY structure is the optimal. On the other hand, qualitative experimental results using a molecular probe (N2 ) suggest that the most stable configuration is OCT [6]. This determination of cluster structure can be, however, sensitive to coverage because cluster-adsorbate interaction may lead to structural changes. Calculations of N2 interaction with small Nin clusters (n ¼ 2; 3; 4) reported a binding energy of 2.1–2.7 eV per N2 adsorption [30]. Therefore, the adsorption of several N2 may release enough energy to produce geometric changes in the Ni6 cluster. Notice that energy differences between isomers and Ni–Ni bond energies shown in Table 1 and Fig. 2, respectively, are smaller than N2 adsorption energy. Another fact that must be taken into account is the possibility that cluster vibrations may cause rearrangements between one nuclear configuration to another, when the temperature increases (a sort of cluster fluxionality). A more detailed analysis of these clusters can be done by evaluating DBE values, bond distances (BD), and charges. Results for PPY, OCT, and CTB are presented in Figs. 2a and c. Values in parentheses correspond to CATIVIC and those in brackets to DFT calculations. BD values of PPY indicate good agreement between CATIVIC and DFT. The distances are slightly distorted in both cases an average value is given for only different types of distances. DBEs
indicate that the most stable bonds correspond to those between the central Ni atom and edge atoms. Charge on the edge atoms are positive, while the central one is negative. In the case of CTP, the predicted CATIVIC structure is more distorted than that obtained by DFT. In both case the shortest distance corresponds to that between the most saturated atoms (Ni(3)–Ni(5)). Note that the longest distance is observed between the less coordinated atoms. The charge on cluster atoms presents the same behavior than in PPY geometry; the most saturated atoms are negatively charged. A distortion is found in OCT structure, particularly CATIVIC results. This is reflected in a small interaction of )7.7 kcal/mol between the shortest vertex distance. Nevertheless, distortion may be explained as a consequence of Jahn-Teller effect that is of relevant importance in geometry structure. A search for different cluster electronic states were not done. However, in some way, MCSA takes into account different electronic states, because there is a correspondence between those and cluster geometries, although the final state depends, in many cases, on the initial-guess density matrix. Values of the HOMO orbital can be associated with the ionization potential. The values for the calculated structures are 6.11, 6.62, and 6.89 eV for OCT, CTP, and PPY, respectively. The correlation with experimental value (6.8 eV) is also a good indication that PPY can be the most stable structure. In addition, recent results from Aguilera et al. [31] suggest, using calculations of average magnetic moment, that Ni6 must have hexahedral geometry, as in PPY geometry. Calculations presented in this work indicate the possibility of using a combination of simple and efficient approaches in obtaining global minimum for metallic clusters where the use of conventional ab initio methodologies becomes prohibitive.
432
F. Ruette, C. Gonzalez / Chemical Physics Letters 359 (2002) 428–433
Fig. 2. Values of bond distances, diatomic bond energies, and charges for PPY, OCT, and CTB isomer Ni6 clusters.
In addition, these results suggest that it is important to search for a global minimum because properties of clusters, such as ionization potential, atomic charges, and magnetism can dramatically depend on the geometry. Our results show the the most stable geometry for Ni6 is not Oh; as it is normally accepted; instead, C5v and C2v symmetries have a lower energy.
Acknowledgements This research has been supported by CONICIT under the contract G-9700667 and Computational Chemistry Laboratory of Physical & Chemical Properties Division of National Institute of Standards Technology for partial supporting the F.R.s sabbatical year.
References [1] J.A. Alonso, Chem. Rev. 100 (1998) 637. [2] M. Ichihashi, T. Hanmura, R.T. Yadav, T. Kondow, J. Phys. Chem. A 104 (2000) 11885, and references therein. [3] N.N. Lathiotakis, A.N. Andriotis, M. Menon, J. Connolly, J. Chem. Phys. 104 (1996) 992, and references therein. [4] J. Guevara, F. Parisi, A.M. Llois, M. Weissmann, Phys. Rev. B 55 (1997) 13283. [5] G.L. Estiu, M.C. Zerner, J. Phys. Chem. 100 (1996) 16874. [6] E.K. Parks, L. Zhu, J. Ho, S.J. Riley, J. Chem. Phys. 100 (1994) 7206. [7] A.N. Andriotis, M. Menon, Phys. Rev. B 57 (1998) 10069. [8] F.A. Reuse, S.N. Khanna, Chem. Phys. Lett. 234 (1995) 77. [9] S.K. Nayak, S.N. Khanna, B.K. Rao, P. Jena, J. Chem. Phys. 101 (1997) 1072. [10] E. Curotto, A. Matro, D.L. Freeman, J.D. Doll, J. Chem. Phys. 108 (1998) 729. [11] M.S. Stave, A.E. Depristo, J. Chem. Phys. 97 (1992) 3386. [12] J.P.K. Doye, D.J. Wales, New J. Chem. 22 (1998) 733.
F. Ruette, C. Gonzalez / Chemical Physics Letters 359 (2002) 428–433 [13] M.A. Nygren, P.E.M. Siegbahn, U. Wahlgren, H. Akeby, J. Phys. Chem. 96 (1992) 3633. [14] M.C. Michelini, R.P. Diez, A.H. Jubert, J. Mol. Struct. 490 (1999) 181. [15] M. Romero, J. Primera, M. Sanchez, A. Sierraalta, S. Brassesco, J. Bravo, F. Ruette, Folia Chim. Theor. Latina 23 (1995) 45. [16] M. Romero, M. Sanchez, A. Sierraalta, L. Rinc on, F. Ruette, J. Chem. Inform. Comp. Sci. 39 (1999) 543. [17] J. Primera, M. S anchez, M. Romero, A. Sierraalta, F. Ruette, J. Mol. Struct. (Theochem) 469 (1999) 177. [18] F. Ruette, M. S anchez, G. Martorell, C. Gonzalez, R. A~ nez, A. Sierraalta, L. Rinc on, C. Mendoza, to be published. [19] G. Blyholder, J. Head, F. Ruette, Theor. Chim. Acta 60 (1982) 429. [20] F. Ruette, C. Gonzalez, A. Octavio, J. Mol. Struct. (Theochem) 537 (2001) 17. [21] J.C. Pinegar, J.D. Langenberg, C.A. Arrington, E.M. Spain, M.D. Morse, J. Chem. Phys. 102 (1995) 666.
433
[22] F. Ruette, F.M. Poveda, A. Sierraalta, J. Rivero, Surf. Sci. 349 (1996) 241. [23] S. Kirpatrick, C.D. Gellat, M.P. Vecchi, Science 220 (1983) 671. [24] D.M. Ferguson, D.G. Garrett, Adv. Chem. Phys. 105 (1999) 311. [25] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [26] P. Roussel-Ragot, G. Dreyfus, in: R. Azencott (Ed.), Simulated Annealing and Parallelization Techniques, Wiley, New York, 1992. [27] C. Gonzalez, in: L.J. Alvarez, (Ed.), Metodos Computacionales en la Catalisis, UNAM, CYTED, Mexico, in press. [28] M.J. Frisch et al., GA U S S I A N 94, Revision D.4, Gaussian Inc., Pittsburgh, PA, 1995. [29] O. Gropen, J. Almlof, J. Chem. Phys. Lett. 191 (1992) 306. [30] F.A. Reuse, S.N. Khanna, B.V. Reddy, J. Buttet, Chem. Phys. Lett. 267 (1997) 258. [31] F. Aguilera-Granja, J.M. Montejano-Carrizales, J.L. Moran-Lopez, Solid State Commun. 107 (1998) 25.