Are zinc clusters really amorphous? A detailed protocol for locating global minimum structures of clusters

Are Zinc Clusters Really Amorphous? A detailed Protocol for Locating Global Minimum Structures of Clusters† Andr´ es Aguadoa,∗, Andr´ es Vegaa, Alexandre Lebonb, and Bernd von Issendorffc Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX First published on the web Xth XXXXXXXXXX 200X DOI: 10.1039/b000000x We report the results of a conjoint experimental/theoretical effort to assess the structures of free-standing zinc clusters with up to 73 atoms. Experiment provides photoemission spectra for Zn− Ncluster anions, to be used as fingerprints in structural assessment, as well as mass spectra for both anion and cation clusters. Theory provides both a detailed description of a novel protocol to locate global minimum structures of clusters in an efficient and reliable way, and its specific application to neutral and charged zinc clusters. Our methodology is based on the well-known hybrid EP-DFT (empirical potential-density functional theory) approach, in which the approximate potential energy surface generated by an empirical Gupta potential is first sampled with unbiased basin hopping simulations, and then a selection of the isomers so identified is re-optimized at a first-principles DFT level. The novelty introduced in our paper is a simple but efficient new recipe to obtain the best possible EP parameters for a given cluster system, with which the first step of the EP-DFT method is to be performed. Our method is able to reproduce experimental measurements at an excellent level for most cluster sizes, implying its ability to locate the true global minimum structures; meanwhile, if exactly the same method is applied based on the existing Gupta potential (fitted to bulk properties), it leads to wrong predicted structures with energies between 1 and 2 eV above the correct ones. Opposite to what was claimed in the past, our work unequivocally demonstrates that Zn clusters are not amorphous, and they rather adopt high symmetry structures for most sizes. We show that Zn clusters have a number of exotic, unprecedented structural and electronic properties which are not expected for clusters of a metallic element, and describe them in detail. 1 Introduction A key step in predicting and rationalising physico-chemical properties of small clusters and nanoparticles is their structural characterization. This is so because the interesting electronic, magnetic, optical, catalytic, etc., properties of nanoparticles depend on the geometry of the ionic skeleton. To achieve a complete characterization would require an exhaustive sampling of the potential energy surface (PES), whose complexity (number of local minima) increases in an exponential way with the number of atoms in the nanoparticle. At low temperature, the most stable structure is that of the global minimum (GM) on the PES, and most experiments probably involve the GM structure, and maybe a few isomers with low excitation energies. As the structure of a free-standing cluster is not di- †Electronic Supplementary Information (ESI) available: Atomic coordinates (in xyz format and ˚ A units) and point group symmetries for the Global Minimum structures reported in this paper; full computational details and benchmark calculations aimed at assessing the accuracy of the level of theory employed. See DOI: 10.1039/b000000x/ a Departamento de F´ ısica Te´ orica, At´ omica y ´ Optica, University of Valladolid, Valladolid 47071, Spain; E-mail: [email protected].uva.es b Laboratoire Chimie ´ Electrochimie Mol´ eculaire et Chimie Analytique, UMR CNRS 6521, 29285 Brest Cedex, France c Physikalisches Institut, Universit¨ at Freiburg, H.-Herder-Str. 3, D-79104 Freiburg, Germany rectly available from experiment, unbiased global optimization (GO) is essential for the prediction and later assessment of cluster structures. The structural assessment itself is usually achieved by indirect means, i.e. by comparing the theoretical predictions on several physical observables with experimental measurements (photoemission and vibrational spectra, as well as diffraction patterns and other properties, are considered as good structural “fingerprints”). Several computational methods for structure optimization of nanoparticles have been devised in recent years. In socalled biased algorithms, the user employs either previously available information on smaller clusters of the same material, or the available structures of chemically similar clusters, or more in general any type of “chemical intuition”, to produce several candidate structures that are then optimized to their respective nearest local minima. These methods can not guarantee the consistent location of the absolute global minimum, as they obviously lack robustness and transferability, although they may be useful in specific cases. Unbiased methods, on the contrary, make no assumptions whatsoever about the structure of the GM, and are the method of choice for a reliable GO search. Most modern unbiased methods integrate an internal local search algorithm as well, which allows the main GO algorithm to operate exclusively with the energies of 1–23 |1 local minima on the PES. Two main types of GO algorithms have been developed. Population-based algorithms deal with a collection of cluster structures, and improve the structure of individuals by moving them on the PES according to some rules that try to mimic biological processes in nature. Iteration of those steps evolves the whole population in such a way that only the best fitted individuals (i.e. those with a lower energy) survive, so that if enough iterations are allowed the probability that the best fitted (GM) individual belongs to the population approaches unity. Evolutionary algorithms, such as the genetic algorithm1,2 and differential evolution3, generate offsprings (new cluster structures) by crossover operations (e.g. “cut-and-splice”) or by mutation of the individuals; the convergence rate of these methods is quite sensitive indeed to the specific set of modifications allowed to the individuals. Particle swarm optimization4,5 rather imitates the movement of individuals in a bird flock, and was first intended for simulating social behavior. The coordinates and velocities of each individual are updated to approach the currently best individual, and the expectation is that the whole swarm will approach in this way the optimal region of the PES around the GM. Individual-based GO algorithms start instead the optimization from a single cluster and allow that individual to sample the whole PES in order to identify the GM structure. As specific examples we mention simulated annealing6, Monte Carlo minimization, minima hopping7,8, stochastic surface walking9,10 and basin hopping (BH)11–14 methods. The essence of the BH method is a mapping of the original PES into a stepwise modified surface obtained by local optimization. In the transformed surface there are plateaus associated with each local minimum, and discontinuous steps separating the different local minima. In this way, energy barriers are substantially decreased so that a Monte-Carlo sampling (the second ingredient in the BH method) is much more efficient than on the original landscape. Several types of movements can be implemented to generate new cluster structures and sample the transformed PES: the basic one is a random change in the atomic coordinates. More specific movements include changes in the coordinates of only surface atoms, or the swap of two different atomic species in the optimization of the chemical ordering of nanoalloys, for example. Although in this paper the focus is on structure optimization using the BH method, we mention in passing that all these methods can also be employed in materials design, where one aims to solve the so-called inverse problem: given a certain desired property, discover (design) a material that possesses that property under a given set of constraints15–18. Most unbiased GO methods require a huge number of evaluations of the energy function in order to achieve an efficient (statistically meaningful) sampling of the PES, and for this reason they have been traditionally tested on simple model surfaces such as those generated by Lennard-Jones, Morse, etc, empirical potentials. But when the goal is to make contact with experimental measurements and deal with realistic systems, sampling an accurate PES becomes indispensable. Then we face an interesting dichotomy and GO methods branch into two different practical implementations. In the first one, the GO method directly samples a quantum accurate description of the PES such as that provided by a density functional theory (DFT) method, for example19,20. This allows to have accurate energies and forces at all stages of the optimization process, but greatly increases the cost of exploring the PES comprehensively. For this reason, the statistical accuracy of the sampling is sacrificed, severely limiting the size of the clusters that it is possible to investigate with confidence. In the second implementation, the GO method is initially applied to an approximate PES generated with a low-cost empirical potential (EP), which allows a statistically meaningful, extensive search. Then a diversity of the structures located in that search are locally re-optimized at the DFT level to recover accuracy in the PES. In this so-called EP-DFT approach21, the potential danger is that one might be re-optimizing at the DFT level structures that have very low stability on the true PES. If we want to obtain results consistent with experiments, it is imperative to train the EP based on first-principles cluster results, a step that is often bypassed simply because there typically are bulk-based parametrizations of the EP available. In this paper, we will show explicitly that transferability of EP parameters is a very serious issue which may lead to unrealistic GM predictions in an EP-DFT search, and will propose a simple and systematic method to fit an EP explicitly to DFT results on metal clusters. We will be concerned specifically with the Gupta potential22–24 but the results should be generally applicable to other EP’s with a relatively small number of parameters. The EP-DFT approach can itself be further subdivided into two quite different practical branches: one may decide to use BH or machine learning-like techniques to train a neuralnetwork potential with a huge number of adjustable parameters in order to ensure high numerical accuracy in the PES generated with the empirical potential25,26. The disadvantage in this method is that the problem becomes a purely numerical one and each parameter lacks a precise physical meaning. We personally prefer to employ an analytic potential with the smallest possible number of parameters compatible with a reasonable level of accuracy in the PES. Even if sacrificing some numerical accuracy, the advantage here is that if few parameters generate reasonable physics, then each parameter has a well-defined physical meaning, which allows to identify useful trends and gain physical insight into the properties of metallic nanoparticles. It is for this reason that we decided to employ the Gupta potential, which has just two independent parameters when expressed in reduced units. As a specific system to test our ideas, we have chosen zinc 2|1–23 clusters. We report putative GM structures for clusters with up to 73 atoms, for both neutral (ZnN) and singly-charged (Zn+ N and Zn− N) clusters, as well as mass and photoemission spectra for cluster anions. One important reason for our choice is that there are many experimental measurements27–29 (including previous mass spectra for cations and our own mass and photoemission spectra for anions) that will serve to validate the theoretical results, and that have not yet been reproduced/interpreted by theory. Photoemission spectra provide very sensitive structural fingerprints, so that a comparison between DFT results for the electron density of states and the photoemission spectra is a recognized method for structure assignment30,31. ZnNclusters are also interesting because in previous theoretical reports they have been predicted to have amorphous GM structures32, which apparently does not match with the well-structured photoemission spectra measured for anions (a well-structured photoelectron spectrum usually indicates an ordered structure, because degeneracies in the electron density of states are related with symmetries of the system). At a fundamental level, ZnNclusters are also interesting as model systems to analyze the insulator-to-metal transition28: while the Zn2dimer and small clusters are weakly bound via vdW effects, the bulk is metallic, so an increasing sp−hybridization and closing of the HOMO-LUMO gap (the gap between the highest occupied and lowest unnoccupied molecular orbitals) is expected as a function of size N. At a more applied level, we are initiating studies on Zn-Mg nanoalloys33 because of their excellent corrosion protection properties. The structures of pure Zn clusters are needed as a reference to calculate excess quantities in Zn-Mg nanoalloys in our future research. 2 Experimental and Computational Methods 2.1 Experiment Zinc cluster anions have been produced in a liquid nitrogen cooled magnetron gas aggregation source, where a pulsed discharge at the end of the aggregation tube was used to enhance the intensity of cluster anions. After exiting the source the clusters were stored in a liquid nitrogen cooled radio frequency octupole trap, where the presence of helium buffer gas leads to an efficient cooling of the clusters. After pulsed extraction from the trap the cluster ions were size-selected in a double reflection time-of-flight mass spectrometer, decelerated and irradiated by 6.4 eV photons (or 4.0 eV in some cases) in the interaction region of a magnetic-bottle-type photoelectron spectrometer. The spectrometer has been calibrated measuring the known photoelectron spectrum of Pt−ions, which leads to an uncertainty in the measured binding energies of less than 30 meV. A subset of the photoemission spectra has been shown in previous reports27,28,41. 2.2 First-Principles calculations The main DFT calculations are performed with the SIESTA code34, under the PBE approximation35 to exchangecorrelation effects. We also employed a van der Waals corrected functional proposed by Klimeˆ set al.36 (KBM functional) in order to check the accuracy of the PBE functional for selected sizes. We include the 3dsemicore states of zinc in the valence active space and describe the effect of the rest of core electrons through a norm-conserving pseudopotential including non-linear partial core corrections. The cluster wavefunction is expanded into a basis of localized atomic orbitals. Benchmark accuracy tests demonstrating the high accuracy of our DFT results, as well as a full description of our computational settings, can be found as Supporting Information (ESI)†. Vertical detachment energies (VDE) of the anions are calculated through a ∆−SCF calculation, that is, as the total energy difference between the anion and the neutral cluster, both clusters adopting the geometry of the anion. In order to simulate the experimental photoelectron spectra, we have broadened each line of the KS eigenvalue spectrum by using a Gaussian function of width 0.10 eV, which leads to a good agreement with experimental line widths. The calculated electronic density of states (EDOS) was then globally shifted in order to align the binding energy of the highest occupied orbital (the HOMO) with the theoretical VDE value. 2.3 Accurate EP-DFT protocol for global optimization Putative GM structures are initially determined through unbiased BH optimizations11 based on a Gupta empirical potential24. We employ the GMIN code developed by David Wales group37,38. For all cluster sizes, we have performed optimizations in which all atoms are given random displacements at each BH step; for some of the bigger sizes, we have additionally performed “seeded” optimizations in which a highsymmetry core is kept fixed and only the rest of atoms is active in the optimization. Around 100 competitive structures are then selected for each size and re-optimized at the KohnSham DFT level. The optimizations for neutral, cation and anion clusters are completely independent runs. In a final “sizecomparison” step, we have always generated new structures for N+1 and N−1 atom clusters by creating a vacancy or adding one adatom to the GM structure of the N−atom cluster. This step is directly performed at the ab initio level. More technical details about the optimization runs can be found in the ESI†. In a combined EP/DFT approach, the accuracy provided by the specific parameterization of the Gupta potential is a crucial issue. No matter how extensive is the basin hopping sampling of the Gupta potential energy surface, nor the intrinsic accuracy of the ab initio DFT model in producing cluster 1–23 |3 energies, the results will be useless (or irrelevant when compared to experimental results) if the DFT energies are calculated from a local optimization of initial structures which significantly depart from the real ones. For this reason, we have carefully tested the accuracy of the current parameterization of the Gupta potential for zinc24 in its application to Zn clusters (we remind the reader that the potential parameters have been fitted only to bulk properties, and transferability is certainly an issue). To simplify the problem, we would like to work with the minimum possible number of relevant parameters, so a first step is simplifying the potential expression as much as possible. The general expression of the Gupta potential energy contains five parameters: E(ZnN) = N ∑ i=1 (Eband i+Erep i)(1) Eband i=−(N ∑ j6=i ξ2exp−2q(ri j r0 −1))1/2 (2) Erep i= N ∑ j6=i Aexp−p(ri j r0 −1)(3) where ri j is the distance between atoms iand j.qand pare dimensionless parameters that determine the effective spatial range of the band and repulsive terms, respectively, while ξ and Afix the strength of those interaction terms and have units of energy. r0is a scaling factor for the distances between atoms, usually identified with the nearest-neighbor distance in the bulk metal. Although the potential incorporates only isotropic interactions, non-additive (many-body) effects are introduced through the square-root dependence of the attractive band energy term on interatomic distances. L´ opez and Jellinek39 demonstrated that only four of the five parameters are truly independent in the Gupta potential expression, and so any one of them, for example r0, can be eliminated in favor of four new independent parameters. Specifically, we define p0=p/r0,q0=q/r0,A0=Aepand ξ0=ξeq. The potential is then rewritten as Ei= N ∑ j6=i A0exp−p0ri j−(N ∑ j6=i ξ02exp−2q0ri j)1/2 .(4) When considering reduced energy and distance units, only two independent parameters remain which determine the true physical content of the Gupta potential40. We have conveniently chosen those parameters as λ=q0/p0=q/pand χ=A0/ξ0= (Aep)/(ξeq). With reduced distances defined as r∗ i j =p0ri j and reduced energies defined as E∗ i=Ei 2ξ0, the reduced total energy of a cluster is E∗=1 2 N ∑ i   N ∑ j6=i χexp[r∗ i j]−"N ∑ j6=i exp[−2λr∗ i j]#1/2   .(5) The two independent parameters λand χcompletely describe the intrinsic topography of the potential energy landscape, including the structure and energetic ordering of local minima (isomers). Later on, in the results section, we will show that these parameters have quite well defined physical meanings, and that their physically meaningful ranges of variation are λ∈(0,0.5)and χ∈(0,∞). Moreover, even if χcan increase without limit from a mathematical point of view, we have explicitly checked that for each fixed value of λ, there is always a maximum value of χbeyond which the topography no longer evolves, essentially corresponding to hard-sphere-like systems. Then, in practice, χ∈(0,χmax(λ)]. We have devised a simple but detailed protocol to choose optimal values of the reduced parameters for zinc clusters, but obviously it could be applied to any other homo-atomic metal clusters. It is based on benchmarking the performance of the reduced Gupta potential against DFT results. There might be several ways to implement this idea in practice, and next we describe ours: First, as we are interested in the size range below 73 atoms, we choose six evenly distributed but otherwise random sizes in the range N=10 −73. DFT calculations aimed to identify the optimal potential will be performed only for those sizes. In second place, and because the parameter space of the reduced Gupta potential is just two-dimensional and bounded, it is a relatively simple and inexpensive matter to perform extensive BH optimizations on a discrete set of (λ,χ)points that provides a representative sample of the whole parameter space. Specifically, we have employed about 100 different potentials with λvalues between 0.01 and 0.49, and with χvalues between 2 and χmax(λ). For each of those potentials, we have performed structure optimizations on the six chosen cluster sizes. The cost of this part of the protocol might be substantially reduced if one is interested in a specific metal for which there is a Gupta potential already available, typically fitted to bulk properties. Then one might explore just a local environment of (λ,χ)values around that reference point which is appropriate for the bulk. Here we have preferred to explore the whole parameter space in order to keep our approach fully unbiased and not to make any assumptions, but we have found a posteriori that the optimal parameters for clusters (see below) are not too distant from the bulk-fitted values. Next, we reoptimize at the DFT level the GM structures located by each potential for the six chosen sizes. A potential is considered more realistic the lower the DFT energy obtained from re-optimization of its GM structure, averaged over the six sizes. Specifically, our objective function is the average binding energy per atom of the six cluster sizes, and the higher this binding energy, the better the potential. This step should clearly show the region of parameter space that is best suited for the considered metal. In our case of Zn clusters, it clearly points to a corner limit of the parameter space, with high val4|1–23 Zn9(D3h)Zn10(C3v)Zn11 + (C3v)Zn11(D3h)Zn11 − (Cs)Zn12 + (Cs) Zn12(Cs)Zn12 − (Cs)Zn13 + (C1)Zn13(C1)Zn13 − (Cs)Zn14 + (Cs)Zn14(Cs) Zn14 − (C1)Zn15 + (C2−C )Zn15(C2v)Zn16 + (C1−C)Zn16 (Cs)Zn16 − (C2v) 2v s Fig. 1 Putative GM structures and approximate point group symmetries of zinc clusters with N=9−16 atoms. The GM structure of charged clusters is shown only when it does not coincide with the GM of the corresponding neutral cluster. When the structure of the N−cluster is obtained by adding one atom to the (N−1)−cluster, the adatom is shown in blue. Blue color is additionally used to highlight that a structure is obtained by adding several atoms to a common structural motif. ues of λand relatively small values of χ. Then we have sampled, in a second step, a finer grid of (λ,χ)points in that local region. In this way, we have found (see next section) that the best agreement with the DFT results for zinc clusters occurs for a set of optimal reduced parameters which, while being located on the same corner of parameter space as the bulkfitted parameters, significantly differ from them, implying the bulk parameters are not transferable to clusters. Once the optimal parameters are located, one can proceed with the usual EP-DFT approach for all sizes. Alternatively, and in order to enhance diversity and so the chances to locate the true GM structure, a good option is to use not only the optimal potential, but a few others in a local environment about the optimal one, to generate candidate structures for DFT re-optimization. Additional technical details about this essential potential optimization step are given in the ESI†. With the new set of parameters, the DFT global minimum is systematically found between the 10-15 most stable Gupta structures for most zinc clusters, which is a very successful result for a potential containing just two parameters. The cluster structures predicted with the bulk parametrization are, on the contrary, between 1 and 2 eV less stable than our putative GM structures when reoptimized at the DFT level, and this is a huge energy error, so the bulk parametrization is simply not representative of the interactions in zinc clusters. This is one of the important new findings reported in this paper, as we have found that any combined EP/DFT approach based on the bulk Gupta parameters is doomed to fail when applied to the location of the putative GM structures of Zn clusters. Accordingly, we include a complete discussion of this issue in the next section. 3 Results and Discussion 3.1 Putative Global Minimum Structures The putative GM structures that emerge from our study are shown in Figures 1-5. A tri-capped trigonal prism (TTP) is obtained for Zn9. The structures in the size range N=10−14 are based on adding atoms to that TTP unit. In some cases, such as Zn+ 11 or Zn+ 12, the adatoms have an anomalously low coordination which is not expected in clusters of a metallic element. This feature, which recurrently appears for some of the bigger size clusters, has been analysed in a recent report,41 where we have shown that those adatoms introduce localized insulating features in an otherwise metallic-like electron density distribution. Therefore, these clusters feature an unprecedented coexistence of two different electronic phases (metallic and insulating) within a single nanoparticle. The GM structure for N=15 is obtained by glueing two TTP units together in such a way that three atoms are common to the two units. The different structures of the cation 1–23 |5 Zn17 + (C2v)Zn17(D5h)Zn17 − (C Zn18 + (D5h)Zn18(C2v)Zn19(C2v)Zn19 − (C3) Zn20(C3)Zn20 − (C3)Zn21(C1)Zn21 − (C1)Zn23(C1)Zn24(Cs)Zn26 + (C1) s) Fig. 2 Putative GM structures and approximate point group symmetries of selected zinc clusters with N=17−26 atoms. Yellow color is used to distinguish the interior or core atoms. Rest of the caption as in Figure 1. and neutral/anion clusters is just due to a different choice of the three common atoms. In fact, the low energy excitations (isomers) identified for Zn15 are just different ways of joining two TTP units, implying that the energy cost of reorienting those two units is much less than the energy needed to destroy one of them. This might lead to interesting thermal behavior at low temperatures, with fluxionality provided by the relative motion of otherwise “rigid” TTP units. The GM structures for size N=16 are essentially obtained by adding one atom to the 15-mer. In the small size range N=9−16, all the atoms still belong to the cluster surface. This is already an indication that zinc clusters are less compactly packed as compared to clusters of other typical metals such as sodium42, aluminum43–45, or silver46, as all those clusters have already developed at least one interior or core atom in this size range. We also notice that the GM structures of neutral, cation and anion are different at many sizes, a feature which will also persist in the bigger clusters discussed below. N=17 is the critical size for the stability of structures without core atoms. Both the cation and the neutral adopt a hollow cage structure, obtained by capping the five square faces of a hollow 12-atom decahedron, and retains full D5hsymmetry in the neutral state. The cation undergoes a slight JahnTeller distortion that lowers the symmetry to C2v. The cage is nearly spherical in shape although it shows a slight prolate distortion, i.e. it is more elongated along its 5-fold symmetry axis than along its equatorial plane. Notice that a “standard” capped decahedron would instead have a significantly oblate shape as the five capping atoms would be more distant from the cage center. However, upon relaxation of that initial structure, the originally square faces of the decahedron elongate along the direction of the 5-fold axis (becoming rectangular facets) while the capping atoms undergo an inwards relaxation and approach the cluster center, resulting in a nearly spherical cage. Hollow cages of a metallic element, that moreover are the GM on the potential energy surface, are very exotic and interesting species as they may be used as hosts for a variety of dopants. In the case of zinc, they are particulary interesting as metallic zinc is widely used in corrosion protecting layers.47 Therefore, zinc nanocoatings may be useful to protect the potentially interesting properties of the embedded material from the environment and other external influences. In a recent report,48 we have shown that Cr@Zn17 is a magnetic superatom with the chromium atom placed in the cage center. This cluster preserves the whole magnetic moment of the isolated chromium atom (6 µB), even after coalescence of two such units or after exposure to an oxygen atmosphere. The GM structure of Zn− 17 is still based on the capped decahedron, but already contains a pseudointernal atom that occupies the center of the decahedron. In the size range N=18 −26 the GM structures contain a single core atom. The 18-mer is obtained by filling the central position of the Zn17 cage. Full D5hsymmetry is preserved only in Zn+ 18, while the neutral and anionic clusters distort and have only C2vsymmetry. An additional capping atom results in the GM structure of Zn19. The GM structure of Zn20 displays a chiral shell with C3symmetry surrounding the internal atom. The stable structures for Zn− 19 and also for sizes N=21 −23 are obtained by removing one atom or adding low coordinated atoms to Zn20. Zn24 features an interesting mixture of structural motifs, based on a decahedron with its equator decorated with trigonal prism units. Zn26 adopts an 6|1–23 Zn28 + (C2) Zn28(C2)Zn28 − (C1) Zn28 2− (D3)Zn29(C2v)Zn29(C2v)Zn29(T) Zn30(C3) Zn30 − (Cs)Zn31 + (C2) Zn31(C2v)Zn32 + (Cs)Zn32(D3h)Zn32 − (C2) Fig. 3 Putative GM structures and approximate point group symmetries of selected zinc clusters with N=28−32 atoms. Rest of the caption as in Figures 1,2. amorphous GM structure. With few exceptions, clusters with N=27 −30 atoms contain an internal dimer. For N=27 −28 the structures have a prolate shape and little symmetry. Zn29 in its neutral state has three nearly degenerate GM structures. One of them is exceptional because it contains a single core atom surrounded by a shell with one of the highest chiral symmetries (T). Cores with two atoms are clearly favored in this size range, and the special stability of the tetrahedral structure is due only to its spherical shape and an electron shell closing effect (Zn29 contains exactly 58 electrons, a magic number of electrons in a spherical jellium picture49,50). The HOMO-LUMO gap of the tetrahedral structure is 1.6 eV, as compared to 0.75 eV for the structures with two core atoms. The near degeneracy is thus the result of a strong competition between electronic and geometric preferences. The average distance between two surface atoms in the T-structure is 2.60 ˚ A, while the average core-shell distance is 3.4 ˚ A. This is just an extreme example of another quite general and defining property of all zinc cluster structures with N≥18, namely the presence of a low-density core, surrounded by a distant, rounded and crowded shell. Because the Zn2− 28 di-anion also contains 58 electrons, we decided to perform unbiased optimizations for its structure and found that it contains also a single core atom (see Fig. 3) and a shell of D3symmetry, resulting from the removal of one atom from the shell of Zn29. It is interesting to notice that all other structures for Zn2− 28 (containing two core atoms) are unstable as the calculated electron affinities of Zn− 28 are negative, and only the D3isomer of Zn− 28 displays a positive electron affinity, which renders Zn2− 28 as the smallest stable zinc di-anion. Its HOMOLUMO gap is 1.2 eV. The two remaining degenerate structures of Zn29 are quite similar to each other. Both contain a decahedral shell with approximate D5hsymmetry that encapsulates a zinc dimer, and differ just in the orientation of that dimer with respect to the cage frame. For Zn+ 29 and Zn− 29, these two structures continue to be nearly degenerate global minima, but the T-isomer is now much less stable. The stability of the T-structure thus depends quite critically on the number of electrons. The GM of Zn30 and Zn+ 30 is obtained by adding one atom to the core of Zn29(T), which lowers the point group symmetry down to C3. This is the second example we find in which the cluster grows by adding atoms to the core rather than to the shell. Zn− 30 already contains three core atoms, although this structure is nearly degenerate with another one obtained by adding a low-coordinated adatom to the C2vstructure of Zn29. The core region of clusters with N=31 −33 atoms contains three atoms in a triangular arrangement, which equips these clusters with an oblate global shape. The way in which the shell grows around this core is highly dependent on the charge state of the cluster and the number of atoms. Zn31, for instance, displays a mixture of hexagonal-closed-packed (hcp) and decahedral growth patterns. Zn32 is a perfect hcp crystalline fragment that retains full D3hsymmetry and so it is the smallest cluster that has the same packing as bulk zinc. In Zn+ 32, the growth of the shell is different on both sides of the triangular core: while on the front side it follows a compact packing pattern, on the back side it adopts a less compact simple hexagonal packing. The size range N=34 −38 features clusters with a tetrahedral 4-atom core. In Zn− 34, the shell preserves the full Td symmetry of the core. This geometric shell closing contains 1–23 |7 Zn33 + (C2v)Zn33(C1)Zn33 − (C1)Zn34 + (C1) Zn36 − (Cs)Zn38 + (C3v)Zn38(C1)Zn38 − (C1)Zn39(C2v)Zn40 + (D2)Zn40(C2v) ) 1 (C 36 Zn) 1 (C 36 Zn+ ) d (TZn34 Fig. 4 Putative GM structures and approximate point group symmetries of selected zinc clusters with N=33−40 atoms. The Tdstructure of Zn34 is shown from two different perspectives. Rest of the caption as in Figures 1,2. four interpenetrating TTP units (the four atoms of the internal tetrahedron being common to the four TTP units). The corresponding neutral cluster adopts the same structure but with a slight distortion that lowers the symmetry down to C2. The HOMO-LUMO gap of Zn34 is only about 0.5 eV, but we observed a large gap between the LUMO and LUMO+1 levels, which motivated us to optimize the Zn2− 34 di-anion. It retains the full Tdsymmetry of the singly-charged anion, and its HOMO-LUMO gap is 1.25 eV, implying that Zn2− 34 is a doubly magic cluster, having closed geometric and electronic shells. Zn+ 38 also displays a highly symmetric (C3v) shell surrounding the tetrahedron, but the rest of GM structures in this size range are of low-symmetry and based on adding atoms to the Tdstructure of Zn− 34. In particular, for Zn35 we obtain a low symmetry structure with a large HOMO-LUMO gap of 1.2 eV, so the electron shell closing for 70 electrons persists even if the geometric shell closing doesn’t. A structural transition in the core region occurs at size N=39 and also for Zn+ 40. Here the core continues to have four atoms, but it is now planar and with a rhombus shape. This transition towards a less compact core allows to pack a larger number of atoms into a single shell of monatomic thickness without modifying the number of core atoms, and it is sometimes observed in zinc clusters as a prelude to the change in the number of core atoms. For both sizes, the shell preserves many of the symmetries of the rhombus core. Clusters with five core atoms are found for most sizes in the interval N=40−44. Zn40 has a multiply twinned layered structure with a trigonal bipyramidal core. Apart from some relaxation of the outermost atoms, this structure can be cut out from a perfect 54-atom decahedron. Zn+ 43 has the same core and a shell that preserves all but the mirror plane symmetries of the core, resulting in a chiral cluster with D3symmetry. Once more, we observed that a transition towards a less compact 5-atom core, in this case a square based pyramid, occurs for example in Zn− 43 and Zn+ 44. Neutral Zn43 and also Zn41 (not shown explicitly) are exceptional in this size range as they feature the reentrance of 4-atom cores. In Zn43, the core is a nearly flat rhombus and the structure contains a hollow site, similar to that observed in Zn17. Six-atom cores occur in clusters with N=45 −48 atoms. Here the reference structure is Zn46, a cluster with 92 valence electrons and so an electronic shell closing according to a spherical jellium model. This cluster is almost perfectly spherical indeed. Its core is octahedral, and the rounded shell maintains a high chiral symmetry (D4) for the whole cluster. Its large HOMO-LUMO gap (1 eV) and its compact shell without defects allows to qualify this cluster as a doubly magic size. All the structures for N=45,47,48 are based on removing or adding atoms to this reference structure. In some cases such as Zn+ 48, we observe again low-coordinated, dangling, adatoms. Clusters with N=49 −53 atoms possess mainly 7-atom octahedral cores of C3vsymmetry, obtained by capping one face of the octahedron, which equips these clusters with an elongated (prolate) shape. Zn49 can still be viewed as grown from Zn46 via addition of three adatoms, but the remaining clusters in this size range do not longer show a resemblance to Zn46. Most of the structures have no symmetries (C1point group) despite there is a clear octahedral order in the core. The only exception is Zn52, for which the shell fully preserves the C3vsymmetry of the core. This structure shows an interesting packing pattern: while the core is octahedral, the 45-atom shell is icosahedral, in fact it is the same 42-atom shell that occurs in a 55-atom icosahedron, with three atoms added on one facet. Mixed packing patterns like this have been found in several metallic nanoalloys,51 but they are not usual in homogeneous metallic clusters. Size N=54 features an 8-atom core with full Tdsymmetry, 8|1–23 43 (D3)Zn+Zn43 (C2v)− 43 (Cs) Zn+ 44 (C2v) Zn46 (D4) Zn+ 47 (C4)Zn+ 48 (C4)Zn49 (C1) Zn51 (C1)Zn52 (C3v)) (C Zn ) Zn 56 Zn (C ) Zn (C ) Zn (C ) Zn65 (C2)Zn+ 66 (C2)Zn− 67 (D3)Zn70 (D3)Zn71 (C2h)Zn73 (C1) 54 (D2d 55 s ) (C2v ZnZnZnZnZn Zn Zn 55 −+ (D3d)56 (D2d)Zn 259 58 s 157 + Fig. 5 Putative GM structures and approximate point group symmetries of selected zinc clusters with N=41−73 atoms. Blue and green colors are used to highlight two different surface facets in Zn70. Rest of the caption as in Figs. 1,2. namely a tetra-capped tetrahedron, surrounded by a complete shell of D2dsymmetry. The high symmetry and the absence of vacancies or adatoms at the shell qualifies this cluster as a geometric shell closing. The GM structures of Zn53 and Zn− 53 (not explicitly shown) are obtained by creating a vacancy on the shell of Zn54, while the 55-atom clusters have one adatom on top of the Zn54 structure. For the Zn− 55 anion, the adatom is located at a higher symmetry site resulting in C2vglobal symmetry and a core with 9 atoms. Zn+ 56 still has an elongated core, namely an octahedron with two atoms capping diametrically opposed facets. This structure has a high D3dsymmetry that is respected by the shell. Neutral and anion clusters with 56 atoms have instead a 10-atom tetrahedral core surrounded by a shell with D2dsymmetry. Zn+ 57 contains one adatom on top of the 56-atom D2dstructure. Different core structures are competitive in this size range, which results in a substantial structural diversity. In the size range N=58 −66, we observe a substantial increase in the rate at which the core size grows. Sizes N=58 −60 have 9-atom cores. For Zn58 the core is octahedral, while for N=59,60 the core is a pentagonal bipyramid with two capping atoms. This core can be viewed as a fragment of a 13-atom icosahedron. Sizes N=61 −63 have no rotational symmetries but contain a 10-atom icosahedral fragment in their cores. For sizes N=64,65 the core is an 11-atom icosahedral fragment. Clusters with 66 atoms and also Zn+ 67 and Zn67 have 12-atom cores. All of them are icosahedral except for Zn+ 66, for which the core is an hcp fragment with stacking sequence 3-6-3. The structures in the core thus seem to be developing a well defined tendency towards icosahedral packing, which is not seen in the free zinc clusters with 9-12 atoms. This is suggestive of important differences in the physical behavior and chemical bonding of zinc atoms in “bulk” (interior) and surface regions. Zn− 67 is the first cluster to possess an icosahedral 13-atom core, and this core is maintained in the whole size range N= 68−71. All these clusters are highly symmetric with the only exception of Zn69 (not shown explicitly), whose shell has no global symmetry. Zn69 has 138 electrons, a magic number in the spherical jellium model. Its shell is indeed very spherical and the cluster has a large HOMO-LUMO gap of about 0.8 eV. Zn− 67 and Zn70, for example, preserve a high chiral symmetry (D3). The growth of the shell onto the icosahedral core is very interesting as it follows different atomic packing schemes on different facets. Some surface facets are Mackay-like (one of these is shown in Fig. 5 in blue color for Zn70); some others adopt a much less compact simple hexagonal packing (green color), so that local TTP units are formed along the 3-fold symmetry axis; the rest of facets adopt an intermediate pattern. In summary, different epitaxies coexist within a single zinc nanoparticle. A similar analysis may be performed for Zn71, which maintains a global C2hsymmetry. Finally, sizes N= 1–23 |9 ric and compact structures with a perfectly icosahedral core). The only exception to this rule is Zn54 which is a geometric shell closing with an open electronic shell, and which is more stable than the electronic shell closure that occurs at Zn56. One reason for this strange behavior is that Zn56 contains 10 core atoms, which is more than optimal for that size. In fact, 10atom cores are not recovered until N=61, so it seems that Zn56 achieves an electronic shell closure at the cost of an unfavourable atomic shell structure. -0,06 -0,03 0 0,03 0,06 Ecoh-Eave (eV) -0,5 0 0,5 ∆2 (eV) -0,5 0 0,5 1 1,5 ∆2 (eV) 10 20 30 40 50 60 70 N -1 -0,5 0 0,5 ∆2 (eV) 10 17 20 15 26 2830 3235,36 25 41 27 46,47 54 54 59 62 65 69 31 34 45 10 1517 20 25 29 32 35 37 41 46 5659 6265 6971 15 19 27 29 31 34 40 46 56 5962 65 68 71 10 20 28 30 35,36 46,47 65 69-71 10 17 20 29 32 35 46 69-71 65 9 15 19 2931 34 46 68-71 39 45 72 54 54 56 54 56 54 Fig. 10 Two different stability measures for Zn clusters. The upper graph shows the PBE cohesive energies, referenced to Eave, which is a four-parameter fit of the form Eave =A0+A1N1/3+A2N2/3+A3N. The three lower graphs show the PBE second-energy differences. Clusters with enhanced stability (or magic sizes) are shown for each of the stability indicators. Black, red and green curves refer to cations, neutrals, and anions, respectively. KBM results are shown for cations in a restricted size range with the dashed black curve, which has been vertically shifted to help visualization. To make connection with experimental abundance spectra, we turn now to an analysis of the second energy differences. First, we compare our results for cations with the mass spectra reported by Diederich et al.29. Marked abundance minima are observed in the experimental mass spectra at N=12 −14,22 −24,29,31,34,37 −39,45,53. All of them are well reflected in our theoretical results either as substantial drops in the cohesive energies or negative values of the second energy difference. The experimentally measured abundance maxima are N=10,(15),20,26,28,30,32,35 − 36,41−42,46−47,54,59,62,65,69,72. The agreement with our theoretical predictions is essentially perfect, suggesting that the identified GM structures are realistic also for cations. A small discrepancy between experiment and PBE results is the maximum predicted by PBE for Zn+ 17, which is absent from the abundance spectrum. To address this issue, we performed systematic KBM calculations (shown with a dashed line) in the restricted size range N=9−21, due to their significantly higher computational expense. Remembering that the KBM and PBE GM structures of Zn− 17 were found to differ (see previous subsection), it is not unexpected to observe some differences between PBE and KBM stabilities. In the case of cations, we find that the KBM and PBE GM structures coincide in this size range, yet the vdW correction to the stabilities is slightly size-dependent, and Zn+ 17 becomes destabilized relative to Zn+ 16 and Zn+ 18. The stability maximum at N=17 dissappears in the KBM results thus improving the agreement with experiment. As compared to the cohesive energies, the second energy differences identify additional magic sizes, that indeed correlate with enhanced experimental abundances for cations. We notice that most of those additional sizes (for example, N=15,41,59,62,65) are located within size intervals of lower-than-average absolute stability, so its abundance in the spectra is only due to their higher local stability, which may have an electronic (N=15,41, see figure 9) or a geometric (N=59,62) origin. We also notice the emergence of apparent “magic number pairs”, which are also a peculiar feature of the experimental mass spectra. These occur for cation sizes bracketing an electron shell closing (which can only occur in neutral clusters). For cations, “magic pairs” are observed for N=10 −11, 20 −21, 32 −33, 35 −36, and 46 −47, as for all those sizes the ∆2value is large and positive. In some cases, the bigger cluster in the pair contains a singlycoordinated adatom with insulating-like features. But in all cases, theory and experiment agree in identifying the smaller cluster as more stable within each pair. In particular, Zn+ 36 and Zn+ 47 stand out as the most stable cations with one unpaired electron outside a closed electronic shell. Abundance mass spectra of anions are shown in Figure 11. Once more, the agreement between theory and experiment is close to perfect: clear maxima in the mass spectrum, many of them followed by sharp dips in the measured abundances, are observed for sizes N= 9,(15),19,(22),25,27,(29),31,34,(40),45 −46,54,56,61 − 62,68,(71). All of them coincide with the highest stabili16 |1–23 112 Chapter 5. Clusters of divalent metals 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 31 25 46 34 27 19 9 3 x10 116 54 68 135 x55 x400 155139 Intensity Cluster size, N Figure 5.9: Mass spectrum of zinc cluster anions. with di↵erent laser energies (h⌫=6.4 eV and 4.02 eV) are presented in the appendix H. The spectra in Fig. 5.10 are cut at a binding energy of 5 eV due to the appearance of two-photon e↵ects (photoelectron spectra of the neutral zinc clusters) at higher BEs for all but the smallest sizes. In most of the spectra one can observe a single peak at low binding energy, which indicates the highest occupied single electron state in the cluster, and a group of peaks at higher binding energy, which represent electron states with more negative total energies. As in the case of magnesium [123] the uppermost peak is assumed to be the lowest p-band state, which is occupied by the single additional electron of the cluster anion, and all other peaks are assumed to be s-band states. The di↵erences between the two highest peaks in the spectra therefore represent the excitation band gaps of the corresponding neutral clusters (in the geometry of the anions) [99]. In Fig. 5.10 strongly structured spectra are observed, suggesting a series of highly degenerate levels (though not necessarily a free-electron shell structure). Most of the cluster sizes up to N= 34 exhibit some band gap at the threshold. This band gap decreases rapidly from size N= 3 to size N= 10; for larger sizes only a slow decrease is visible. Figure 5.11 shows the size dependence of the band gap. The existence of a significant gap at least up to size N=25 (51 electrons) indicates that zinc clusters in this size range have not fully undergone the transition to a metal and therefore behave like semiconductors. Fig. 11 Abundance mass spectrum of cluster anions, Zn− N. ties predicted by theory. Magic pairs are now expected at N=9−10,19 −20,31 −32,34 −35,45 −46,68 −69, and many of them are indeed observed both in theory and experiment. In most cases we observe again that the smaller cluster in the pair is clearly more stable, and moreover this assymetry in the stabilities around a closed shell is more pronounced than in cations. Zn− 46 is the only exception to the general rule, as it is more stable than Zn− 45. This points to a substantial contribution of the geometric shell closing to the stability of that particular cluster size, as it is the only magic number which is the same in the cation, neutral, and anion stabilities. This is a sensible result as Zn46 is a doubly magic cluster, having simultaneously closed electronic and geometric shells. A similar double-closure effect occurs for Zn29 and Zn32 in the neutral series and would be expected for Zn2− 34 in the di-anion series. 4 Performance of the Gupta potential A final interesting point in our discussion is that the gross structural features of most clusters are already captured by the Gupta potential. This is of course rewarding, but at the same time is strange: as explained in our previous report,41 the chemical bonding at the shell of zinc clusters displays a substantial directionality which is reminiscent of an sp2−hybridization mechanism, and the Gupta potential certainly does not include any anisotropic interactions, at least not in an explicit way. It is possible, however, that the many-body character of the interactions can produce some effective directionality effects. As shown in the theory section, the energy expression depends exclusively on two independent parameters when expressed in reduced units, and we conveniently choose those parameters as λ=q/pand χ=Aexp(p) ξexp(q). We would like to analyze the physical meaning of those parameters in order to explain why is the Gupta potential able to predict these exotic structures for metallic clusters. 510 15 Reduced Distance -1 0 1 Reduced Energy λ=0.1 λ=0.2 λ=0.3 λ=0.4 λ=0.5 0510 15 Reduced Distance χ=100 χ=50 χ=20 χ=10 χ=5 χ=1 0510 15 20 Reduced Distance Na Cu Al Au Pb Zn(b) Zn(c) χ=100 λ=0.4 Fig. 12 The left and middle graphs show, respectively, the λ−and χ−dependences of the reduced Gupta energies of a perfectly square tetramer as a function of the nearest-neighbor reduced distance. The right graph shows similar plots for selected metals24,60. Zn(b) and Zn(c) refer to the bulk and cluster fitted potential for zinc, respectively. As a simple test system, we have chosen a perfectly square tetramer which already includes many-body interactions but can be described in terms of a single nearest-neighbor distance. We show in Figure 12 the total reduced energy of such a tetramer as a function of the nearest-neighbor reduced distance, for a variety of parameter values. On the left graph we show the energies obtained with a common value of χ=100 and different values of λ. These energy curves display very similar repulsive cores at short distances, and differ most markedly in the attractive and tail regions of the potential. The parameter λ=q/p, which by definition measures the spatial range of the attractive (many-body) part of the energy relative to the range of the pairwise repulsion, essentially controls the tail region of the potential. Therefore, with χfixed, λwill mainly determine the packing density or average coordination number: lower λvalues imply significant attractive interactions between more distant atoms, which will tend to produce strained structures with contracted radial distances and expanded tangential distances, a large internal radial pressure, a high average coordination number, and a low ratio of surface to core atoms. These features are typical of icosahedral or poly-tetrahedral structures, for instance. As λincreases the many-body term becomes of shorter range, and GM structures develop a lower average coordination number, i.e. they become less densely packed. The large radial pres1–23 |17 sure continuously decreases so less strained structures such as decahedra, and later non-strained structures, emerge. For very large λ−values the optimal packing density around internal atoms is low enough that a significant inverted strain can be established, leading to expanded radial distances and contracted tangential distances, with a high ratio of surface to core atoms. The middle plot shows the energies obtained with λ=0.4 and different values of χ. The specific value of λcorresponds to a rather short-range potential but the general conclusions extracted concerning the physical meaning of χdo not depend on the value chosen for λ. Now the different curves tend to the same long-distance energy and differ most significantly at short distances. χis defined as the ratio of the repulsive to the attractive energies at ri j =0, and the plot demonstrates that it essentially controls the core part of the potential. At fixed λ, higher values of χincrease the stickiness of the interactions and thus reduce the ability of the cluster to tolerate bond strain. Lower values of χproduce softer-core potentials which admit an increasingly larger amount of strain. For any fixed value of λ, we therefore expect a transition from strained structures (at low χ) to non-strained structures such as fcc fragments as the value of χincreases, the specific λvalue mostly determining the packing density, i.e. the ratio of surface to core atoms. The intervals of physically meaningful values for the two parameters are χ>0 (obvious; otherwise there is no repulsion to balance the attractive force) and λ∈(0,1/2). The reason for the restriction λ<1/2 or equivalently q<p/2 is not so obvious. In a simplified lattice model that assumes no strain (all atoms at the same interatomic distance) and neglects nonnearest-neighbor interactions, it can be easily shown61 that the dependence of the Gupta total energy on the average coordination number nis dE dn =p−2q 2(p−q) E n.(6) At the critical value λ=1/2, the energy becomes independent of the coordination number, and for λ>1/2 the potential favors the lowest coordination, so the optimization directly leads to fragmentation of any cluster into dimers. Although this calculation corresponds to a simplified model, the true system approaches indeed that simplified limit as λapproaches its critical upper value: the interactions become of such a short-range that non-nearest-neighbor interactions are negligible, and the strain is so strongly penalized that the assumption of no strain is essentially correct. The two parameters λand χcan be independently fixed in the Gupta potential, which provides additional flexibility over, for instance, a Morse potential that contains a single free parameter when expressed in reduced units62. Thus, in the Gupta potential the softness of the core region and the spatial range of the attractive part are independent variables. The interplay between them provides new insights on our present understanding of molecular potential energy landscapes. For example, the shorter the range of a Morse potential the larger the number of local minima, featuring a potential energy landscape that is globally smooth but locally rough, an effect that can be explained using catastrophe theory.62 On the contrary, in a Gupta potential of arbitrarily short range (for example, λ=0.49), we have observed that the number of local minima decreases upon softening the core repulsion, so it is always possible to generate locally smooth energy landscapes with few local minima even for short-range potentials. An interesting observation is that the critical χmax value beyond which only non-strained structures are stable becomes larger the smaller the value of λ. For example, for λ=0.05 the whole range between χ=1 and χmax =1080 is physically meaningful, while for λ=0.49 the interval of physically different potentials narrows down essentially to χ∈(0,100]. In particular, for very large λ−values, an apparently small change in χcan produce a completely different potential and, of course, very different GM structures; for low λ−values, a much larger change of χis needed in order to observe significant differences in the predicted structures. The rightmost graph in Figure 12 shows the reduced Gupta energies of a variety of metals for which a parameterization already exists in the literature.24,60 Different metals are located in different regions of the two-dimensional parameter space defined by the (λ,χ) pair. The potential for sodium, for example, incorporates a long-range attraction (λ=0.128), and with χ=374.7 is not sticky at all. In this range of parameter space, the GM structures include strained poly-tetrahedral and icosahedral motifs with crowded cores, shorter radial distances and longer distances between surface atoms42. The aluminum potential (λ=0.292, χ=41.20) is of shorter range but still tolerates strain quite well, resulting in a competition between icosahedral and decahedral structures45. Gold (λ=0.39, χ=56,3) and lead (λ=0.38, χ=40.2) both display sticky and short-range interactions, which can produce a tendency towards amorphization.61 The parameters for zinc, fitted to bulk properties24, are λ=0.46, χ=26.87, and this potential also tends to produce amorphous cluster structures.32 In this work, we have found that λ=0.48, χ=2 are (on average) the most appropriate parameters for zinc clusters, i.e. those that produce the GM structures that best agree with the DFT results. The modification of the χ−parameter as compared to the bulkfitted value is substantial, as the new potential features a much softer core. We notice that the cluster structures previously obtained with the bulk parametrization32,53 are, according to our DFT calculations, between 1 and 2 eV less stable than the GM structures predicted in this paper. Moreover, those structures were predicted to be amorphous while the improved Gupta potential results in high-symmetry structures for many sizes in agreement with the DFT results, emphasizing the importance 18 |1–23 of fine-tuning potentials for cluster studies. And most importantly, only the present results can match the photoemission measurements and the abundances observed in mass spectra. The new zinc potential belongs to a region of the parameter space that apparently has not been previously explored. It has an extremely short range, the shortest one of all metals shown in the figure. But at the same time, the potential has an extremely soft core. The combination of both factors produces strained structures with crowded and rounded shells and lowdensity cores, that agree very well with the DFT results. These structures have longer radial distances and shorter distances between surface atoms, just the opposite trend as compared to the more typical and better characterized poly-tetrahedral, icosahedral or decahedral structures. Strain is allowed by the very soft core, even if this is unusual in short-range potentials. On the other hand, the extremely short-range nature of the attractive interactions favors a low packing density around core atoms and thus pushes the shell outwards. Alternatively, and using the jargon of embedded atom model (EAM) potentials, a short range attraction implies that atoms achieve their optimal embedding density already at low coordination. In summary, through its coordination preferences, the isotropic Gupta potential favors a long separation between core and shell and thus effectively induces a tendency towards twodimensionality. The sp2−like hybridization observed within the shell at the DFT level is just the response of the electron density (only available from the DFT results) to that specific ionic potential. We close this section by explicitly noticing that the Gupta potential does not predict singly-coordinated adatoms for any size. This is what should be expected, because it is designed to describe metallic bonding and we have demonstrated that adatoms involve the local re-entrance of insulating bonding features. Therefore, the pool of trial structures generated by the BH sampling for clusters such as Zn21, for example, may not be considered a good starting point for the DFT reoptimizations. Fortunately, we have found that all those clusters are extremely good structure seekers. For Zn21, even if none of the trial structures contained an adatom, in more than 90% of the local DFT optimizations an atom was expelled from the cluster and promoted to an outer shell as an adatom. This observation implies that for those specific sizes, all the structural isomers without adatoms represent very high energy excitations (typically more than 1 eV less stable than the DFT GM). Essentially, no matter how good is the starting trial structure, it almost always converges to an adatom structure. The first-principles calculations are completely necessary to describe the reentrance of insulating behavior after the electron shell closings. 5 Conclusions and Summary We have performed a conjoint theoretical/experimental study aimed at determining the GM structures of neutral and singlycharged zinc clusters with up to 73 atoms. Experiment provides photoemission spectra as structural fingerprints for cluster anions, as well as mass spectra for both anion and cation clusters. Theoretical calculations are based on an EP-DFT global optimization approach, in which the approximate energy landscape generated by a Gupta potential model is first sampled with basin hopping optimizations, and then several (around 100) of the more stable isomers so identified are locally reoptimized on the first-principles energy landscape predicted by DFT. The EP-DFT approach is both convenient and necessary for the relatively large clusters studied in this paper, as the statistical complexity of the potential energy surface makes impossible a direct sampling at the DFT level. A key point emphasized in this manuscript is that the success of the EP-DFT approach critically depends on the ability of the EP to reproduce, at a reasonable level at least, the essential topographical features of the DFT potential energy surface. Therefore, training the EP to reproduce DFT properties of clusters should be in our opinion an inescapable step in any EP-DFT approach, even though this step is by-passed in most previous papers employing this method (which directly employ EP parameters fitted to bulk properties). In this paper, we have proposed an original detailed protocol for deducing the optimal EP parameters for cluster studies. First, a redundant parameter of the Gupta potential is eliminated and then the potential is expressed in reduced units, where it has only two independent parameters. Our proposal is then simple: the optimal potential in the two-dimensional parameter space is that one which produces GM structures which are most stable according to DFT. As simplistic as it may look like, EP-DFT optimizations carried out with this proposal are able to identify the correct GM structures of zinc clusters for nearly all the sizes considered, as demonstrated by the excellent reproduction of the photoemission spectra for anions and also of the mass spectra for both cations and anions. Exactly the same EP-DFT optimizations, but carried out with an EP fitted to bulk properties, produces instead wrong GM structures which are between 1 and 2 eV less stable than the ones identified here. Failing to train the EP specifically for clusters may therefore produce huge quantitative and qualitative errors. In the case of zinc clusters, the previous structural EP-DFT prediction, performed with an EP fitted to bulk properties, produced amorphous structures,32 incompatible with the highly structured photoemission spectra for most sizes. Our work clearly demonstrates that only a few clusters of zinc are amorphous; most of them have instead a high structural order and so highorder point group symmetries. Our protocol can be employed 1–23 |19 for predicting the GM structures of other clusters (metallic or not) and in conjunction with other empirical potentials, so it has quite general implications for any future work focused on structure prediction. It is, for example, quite probable that previously reported EP-DFT structures for other metal clusters, obtained with a bulk-fitted and non-transferable EP, carry substantial errors as well, if there are no experimental fingerprints available for a definitive structural assessment. Our method thus opens the way to systematic revisions of medium-sized metallic cluster structures in future studies. Another interesting conclusion from our work is that a simple potential with only two independent parameters in reduced units can describe very well the exotic ordered structures of zinc clusters, only if properly trained in the way described in this paper. The lesson to be learnt is that taking this step can really be worth the effort, and should be done before blaming a given potential and before trying a more complex potential form with more independent parameters. Our approach pushes the Gupta model to its best possible performance level for a given cluster system. Expressing the Gupta potential in reduced units before training it is an essential step in our approach that must be highlighted. It is not only convenient from a computational point of view as it obviously reduces the number of parameters to be fitted; it rather has additional physical implications. Using the relations r∗ i j =p0ri j and E∗ i=Ei 2ξ0, any reduced Gupta potential can be re-scaled to obtain a potential with full units and the desired average distance and cohesive energy. We have explicitly checked that the Gupta potential is unable to describe the correct size dependences of distances and binding energies in the small non-scalable size regime considered in this paper. In words, the p0and ξ0values needed to recover the right average distance and binding energy of, say, Zn20, are different from those obtained for Zn40, so if working with full units, a different Gupta potential would be needed for each size, while in reduced units a single potential is uniformly good in the whole size range below ∼100 atoms!. Trying to fit the full-units potential on DFT data for a variety of cluster sizes would amount to request from the potential something it can not describe, and it would introduce significant noise and innacuracy in the fitting. Moreover, absolute distances and energies are not relevant for the EP part of an EP-DFT approach, as the DFT calculation will fix the correct values of them in the end. Fitting the reduced potential in the way proposed here has the virtue of focusing on the essential topographical features of the energy landscape, including the hierarchy of local minima energies and avoiding the noise coming from the absolute scales of the potential energy surface. We believe this is the key to the success of our approach. An interesting issue which we leave for future work is the performance of our potential in the bulk limit (at the moment our code does not contemplate periodic boundary conditions and it will take additional work for it to deal with bulk systems). The parameters reported by Cleri and Rosato24 were fitted to bulk values of lattice constant, cohesive energy and elastic constants for each metallic element. As we have seen, the lattice constant and cohesive energy contribute no information at all about the topography of the energy surface, and only describe the less interesting absolute scales, i.e. all reduced potentials might be re-scaled to right values of those two quantities. So only the elastic constants remain to contribute some information about the topography, and that information is very local in nature, coming from small displacements around a single minimum. We have explicitly checked that our potential, when re-scaled to full units, reproduces accurately the DFT frequency of a monopolar “breathing” mode for Zn clusters, suggesting it may provide reasonable vibrational frequencies. We believe that our potential may provide less accurate elastic constants but more accurate relative stabilities of different crystalline phases, though this remains at the moment as a pure speculation. We have applied our EP-DFT protocol to describe the main structural properties of zinc clusters. Our results show that most zinc clusters are not amorphous, and rather adopt highsymmetry structures that exhibit some exotic properties for a metal. First, the structures are less compactly packed as compared to other typical metals, meaning that the number of interior atoms is quite lower. The most extreme example is Zn17, whose structure is a hollow cage. Interestingly, similar core-shell structures with an anomalously small number of core atoms have been suggested for some medium-sized gold clusters as well63. This cage can host dopant atoms at its center while having little interaction with them, and in the case of transition metal dopants the cage is able to maintain a high magnetic moment.48 Second, the clusters display an anomalously long radial separation between the few atoms in the core and a crowded rounded shell of atoms. This means that zinc clusters have an exotic strain distribution as compared to most other metallic clusters, displaying contracted tangential distances at the crowded shell and expanded radial distances between core and shell atoms. This is exactly opposite to the strain distribution in icosahedral or poly-tetrahedral clusters, for example. The long core-shell distance has interesting structural and dynamical consequences. At the structural level, it implies that there are no strong epitaxial preferences for the growth of the shell onto the core, and in fact we have found a coexistence of several epitaxial growths in several clusters, as for example Zn70. At the dynamical level, it implies that there are low energy barriers opposing the relative rotation of the core with respect to the shell. For example, the dimer core in Zn29 can rotate almost freely inside the decahedral cage,41 and in general core-shell relative rotations always appear between the lowest-energy vibrational excitations. If the Zn2dimer could be substituted with some 20 |1–23 magnetic transition metal dimer, and if the property of nearlyfree rotation of that dimer could be preserved, it would generate a “nano-compass”, a possibility we are presently exploring. Third, many of the highly-ordered GM structures are chiral, the most symmetric one being Zn34 with tetrahedral chiral symmetry, which suggests future applications of Zn clusters in optical activity studies. Fourth, there are some exceptional cluster structures following an electron shell closing, which contain a singly-coordinated adatom (examples are found at N=11,12,21 and 48 and are more prominent in cluster cations). Another exceptional cluster though for a different reason is Zn32, which is the first and only cluster that is a fragment of the bulk hcp lattice. Finally, we mention that the GM structure is strongly dependent on the number of electrons, as for many sizes the structures of neutral, cation and anion differ. Dispersion or van der Waals effects have been found to have a small but systematic effect, namely tending to increase the number of core atoms. The structures reproduce very well the experimental measurements for nearly all sizes. We have also described the electronic structure of zinc clusters. Most clusters are very well described by the jellium model predictions, and can thus be viewed under a superatom metallic picture. This can be surprising at first sight given the non-metallic structural features identified in the previous paragraph. First, the long core-shell separation produces a deep minimum of electron density along the radial direction, and in a previous work41 we have shown that this minimum is as deep as bewteen two ions of a typical insulating material such as LiF. This is admittedly not a jellium feature, but in practice we have demonstrated here that it provides only a weak perturbation on top of the jellium picture. Specifically, high angular momentum levels are stabilized relative to low angular momentum levels because they better overlap with the crowded ionic density in the shell, but the stabilization effect is small and does not produce unexpected magic numbers or substantial reordering of jellium levels. The only exception is hollow Zn17 for which the 1Fshell is more stable than the 2Sshell. The conclusion to extract is that angular delocalization is the essential property for the validity of the superatom picture, while radial delocalization is not so important. Most of the magic numbers can thus be explained by jellium predictions, and the rest of magic numbers could be related to geometrical shell closings. Second, singly-coordinated adatoms are also not an expected metallic property. In a previous work41 we have shown that in those clusters an insulating bond between adatom and superatoms coexists with the metallic bonding within the superatom. Accordingly, the two localized electrons in the insulating bond populate an “impurity-like” level within a jellium density of electronic states, and the jellium picture continues to be valid to a large extent. The consequence is that clusters with two electrons more than a jellium shell closing continue to have a large HOMO-LUMO gap and so are also shell closings, with an impurity level populated at an energy quite below the HOMO. The size evolution of zinc clusters towards full metallicity is therefore quite gradual even if most clusters display typically metallic features. We gratefully acknowledge the support of the “Junta de Castilla y Le´ on” (Project No. VA124G18). Facilities provided by the Pole de Calcul Intensif pour la Mer, DATARMOR (Brest) are also acknowledged. References 1 Deaven, D. M.; Ho, K. M. Phys. Rev. Lett. 1995,75, 288291. 2 Johnston, R. L. Dalton Trans. 2003, 4193-4207. 3 Storn, R,; Price, K. J. Global Optim. 1997, 341-359. 4 Zhang, J.; Dolg, M. Phys. Chem. Chem. Phys. 2015,17, 24173-24181. 5 Avenda˜ no-Franco, G.; Romero, A. H. J. Chem. Theory Comput. 2016,12, 3416-3428. 6 Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Science 1983, 220, 671. 7 Goedecker, S. J. Chem. Phys. 2004,120, 9911. 8 Amsler, M.; Goedecker, S. J. Chem. Phys. 2010,133, 224104. 9 Shang, C.; Liu, Z. P. J. Chem. Theory Comput. 2013,9, 1838-1845. 10 Zhang, X. J.; Shang, C.; Liu, Z. P. J. Chem. Theory Comput. 2013,9, 5745-5753. 11 Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997,101, 5111. 12 Wales, D. J.; Scheraga, H. A. Science 1999,285, 13681372. 13 Oakley, M. T.; Johnston, R. L.; Wales, D. J. Phys. Chem. Chem. Phys. 2013,15, 3965-3976. 14 Calvo, F.; Schebarchov, D.; Wales, D. J. J. Chem. Theory Comput. 2016,12, 902-909. 15 Dieterich, J.; Hartke, B. Phys. Chem. Chem. Phys. 2015, 17, 11958-11961. 16 Cerqueira, T. F. T.; Sarmiento-P´ erez, R.; Amsler, M.; Nogueira, F.; Botti, S.; Marques, M. A. L. J. Chem. Theory Comput. 2015,11, 3955-3960. 17 Yu, M.; Yankovich, A. B.; Kaczmarowski, A.; Morgan, D.; Voyles, P. M. ACS Nano 2016,10, 4031-4038. 18 Fournier, R.; Mohareb, A. J. Chem. Phys. 2016,144, 024114. 19 Davis, J. B. A.; Shayeghi, A.; Horswell, S. L.; Johnston, R. L. Nanoscale 2015,7, 14032-14038. 20 Shayeghi, A.; G¨ otz, D.; Davis, J. B. A.; Sch¨ afer, R.; Johnston, R. L. Phys. Chem. Chem. Phys. 2015,17, 2104-2112. 1–23 |21 21 Ferrando, R.; Fortunelli, A.; Johnston, R. L. Phys. Chem. Chem. Phys. 2008,10, 640-649. 22 R. P. Gupta, Phys. Rev. B 1981, 23, 6265. 23 V. Rosato, M. Guillop´ e, B. Legrand, Philos. Mag. A 1989, 59, 321. 24 Cleri, F.; Rosato, V. Phys. Rev. B 1993,48, 22. 25 Ouyang, R.; Xie, Y.; Jiang, D. Nanoscale 2015,7, 1481714821. 26 Behler, J. J. Chem. Phys. 2016,145, 170901. 27 Kostko, O.; Wrigge, G.; Cheshnovsky, O.; von Issendorff, B. J. Chem. Phys. 2005,123, 221102. 28 von Issendorff, B.; Cheshnovsky, O. Annu. Rev. Phys. Chem. 2005,56, 549-580. 29 Diederich, Th.; D¨ oppner, T.; Fennel, Th.; Tiggesb¨ aumker, J.; Meiwes-Broer, K.-H. Phys. Rev. A 2005,72, 023203. 30 Binggeli, N.; Chelikowsky, J. R. Phys. Rev. Lett. 1995,75, 493. 31 Akola, J.; Manninen, M.; H¨ akkinen, H.; Landman, U.; Li, X.; Wang, L.-S. Phys. Rev. B 2000,62, 13216. 32 Michaelian, K.; Beltr´ an, M. R.; Garz´ on, I. L. Phys. Rev. B. 2002,65, 041403. 33 Lebon, A.; Aguado, A.; Vega, A. Corros. Sci. 2017,124, 35-45. 34 Soler, J. M.; Artacho, E.; Gale, J. D.; Garc´ ıa, A.; Junquera, J.; Ordej´ on, P.; S´ anchez-Portal, D. J. Phys.: Condens. Matter 2002,14, 2475. 35 Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996,77, 3865. 36 Klimeˆ s, J.; Bowler, D. R.; Michaelides, A. J. Phys.: Condens. Matter 2010,22, 022201. 37 http://www-wales.ch.cam.ac.uk/GMIN/. 38 D. J. Wales, “Energy Landscapes and Structure Prediction Using Basin-Hopping” in “Modern Methods of Crystal Structure Prediction”, edited by A. R. Oganov, WileyVCH, Weinheim, 29–54, 2010. 39 L´ opez, M. J.; Jellinek, J. J. Chem. Phys. 1999,110, 8899. 40 We notice that reduced units were already employed in the past, but without ellimination of redundant parameters, in Garz´ on, I. L.; Jellinek, J. Z. Phys. D 1991,20, 235-238; Z. Phys. D 1993,26, 316-318. 41 Aguado, A.; Vega, A.; Lebon, A.; von Issendorff, B. Angew. Chem. Int. Ed. 2015,54, 2111-2115. 42 Aguado, A.; Kostko, O. J. Chem. Phys. 2011,134, 164304. 43 Starace, A. K.; Neal, C. M.; Cao, B.; Jarrold, M. F.; Aguado, A.; L´ opez, J. M. J. Chem. Phys. 2009,131, 044307. 44 Ma, L.; Issendorff, B.; Aguado, A. J. Chem. Phys. 2010, 132, 104303. 45 Aguado, A.; Jarrold, M. F. Annu. Rev. Phys. Chem. 2011, 62, 151. 46 Chen, M.; Dyer, J. E.; Li, K.; Dixon, D. A. J. Phys. Chem. A2013,117, 8298. 47 Leygraf, C. Atmospheric Corrosion,2000, Wiley, New York. 48 Lebon, A.; Aguado, A.; Vega, A. Phys. Chem. Chem. Phys. 2015,17, 28033-43. 49 Brack, M. Rev. Mod. Phys. 1993,65, 677. 50 de Heer, W. A. Rev. Mod. Phys. 1993,65, 611. 51 Aguado, A. in “Nanoalloys: From Fundamentals to Emergent Applications”, edited by Florent Calvo, 2013, Elsevier. 52 Aguado, A. J. Phys. Chem. C 2012,116, 6841-6851. 53 Doye, J. P. K. Phys. Rev. B. 2003,68, 195418. 54 Wang, J.; Wang, G.; Zhao, J. Phys. Rev. A. 2003,68, 013201. 55 Li, B.; Zhu, Y.; Lu, C.; Ye, G. J. Clust. Sci. 2017,28, 3281-3298. 56 Iokibe, K.; Tachikawa, H.; Azumi, K. J. Phys. B: At. Mol. Opt. Phys. 2007,40, 427-436. 57 Oymak, H.; Erkoc, S. Int. J. Mod. Phys. B 2012,26, 1230003. 58 Gutsev, G. L.; Weatherford, C. W.; Belay, K. G.; Ramachandran, B. R.; Jena, P. J. Chem. Phys. 2013,138, 164303. 59 Becke, A. D.; Edgecombe, K. E. J. Chem. Phys. 1990,92, 5397-5403. 60 Li, Y.; Blaisten-Barojas, E.; Papaconstantopoulos, D. A. Phys. Rev. B 1998,57, 15519. 61 Soler, J. M.; Beltr´ an, M. R.; Michaelian, K.; Garz´ on, I. L.; Ordej´ on, P.; S´ anchez-Portal, D; Artacho, E. Phys. Rev. B 2000,61, 5771-5780. 62 Doye, J. P. K.; Wales, D. J. J. Chem. Soc. Faraday Transactions 1997,93, 4233-4243. 63 Shao, N.; Huang, W.; Gao, Y.; Wang, L.-M.; Li, X.; Wang, L.-S.; Zeng, X. C. J. Am. Chem. Soc. 2010,132, 65966605. 22 |1–23 Zn34 Td Zn46 D4 Zn54 D2d −− − Fig. 13 TOC figure. TOC text: “A novel computational protocol demonstrates that Zn clusters are not amorphous and reproduces photoemission and mass spectra.” 1–23 |23