Study on Structural Evolution, Thermochemistry and Electron Affinity of Neutral, Mono- and Di-Anionic Zirconium-Doped Silicon Clusters ZrSin0/-/2- (n = 6–16)

We have carried out a global search of systematic isomers for the lowest energy of neutral and Zintl anionic Zr-doped Si clusters ZrSin0/-/2- (n = 6–16) by employing the ABCluster global search method combined with the mPW2PLYP double-hybrid density functional. In terms of the evaluated energies, adiabatic electron affinities, vertical detachment energies, and agreement between simulated and experimental photoelectron spectroscopy, the true global minimal structures are confirmed. The results reveal that structural evolution patterns for neutral ZrSin clusters prefer the attaching type (n = 6–9) to the half-cage motif (n = 10–13), and finally to a Zr-encapsulated configuration with a Zr atom centered in a Si cage (n = 14–16). For Zintl mono- and di-anionic ZrSin-/2-, their growth patterns adopt the attaching configuration (n = 6–11) to encapsulated shape (n = 12–16). The further analyses of stability and chemical bonding make it known that two extra electrons not only perfect the structure of ZrSi15 but also improve its chemical and thermodynamic stability, making it the most suitable building block for novel multi-functional nanomaterials.


Introduction
Transition metal (TM) silicides as novel functional materials have been extensively utilized in diverse fields such as aerospace, microelectronic device manufacturing, optical instruments, magnetic material, and the chemical industry because of their oxidation resistance, high-temperature stability, and low electrical resistivity [1,2]. For example, in the early 1960s platinum silicide was successfully used to improve the rectifying characteristics of diodes [3]. Nickel silicide showed significantly enhanced catalytic activity in CO methanation [4]. Composite coatings including molybdenum silicide can be used for protection of missile nozzles and aero engine blades [5]. Zirconium silicide is an ideal device with significant potential applications as anode materials in Li-ion batteries, protective coating on zirconium-alloy fuel cladding in light water reactors, heating elements for electrical furnaces, and also can be incorporated into the semiconductor industry as large scale integrated circuits, Schottky barriers, Ohmic contacts, and rectifying contacts due to the fact that interfaces of zirconium and silicon have good lattice matching, sharp interfaces, low Schottky barrier heights, high conductivity, and excellent thermal stability [1,2,[6][7][8]. Therefore, a large amount of research on zirconium silicides investigations of neutral and anionic Zr-doped silicon clusters are deficient, especially for the important Zintl anions ZrSi n -/2-. As a matter of fact, only by understanding the ground state structure and the growth pattern of Zr-doped Si clusters can we explain the diversity of their properties. Herein, the unbiased ABCluster global search technique coupled with the mPW2PLYP method and the cc-pVTZ basis set for Si atoms and the cc-pVDZ-PP basis set for Zr atoms is utilized for the configuration optimization of neutral and Zintl anionic ZrSi n 0/-/2-(n = 6-16) with the objective of systematically investigating their structural stability and growth patterns, accurately probing the thermochemistry and electron affinity, illuminating details of the bonding characteristic, and offering significant data for further theoretical and experimental investigations of silicon clusters doped with other TM atoms.

Evolution of the Ground State Structure
The global minimal structures and/or competing global minimal isomers optimized with the mPW2PLYP/MIDDLE method for neutral, mono-, and di-anionic ZrSi n 0/-/2species are illustrated in  Table S1.
The average binding energy (E b ), electronic state, and NPA charge on Zr atoms for the global minimum are provided in Table 1. The ground state for all mono-anionic ZrSin -(n = 6-16) is calculated to be a doublet. Their structural growth pattern prefers the attaching type to the encapsulated shape with the Zr atom centered in the Si cage. The global minimum of ZrSi6 -is a distorted pentagonal bipyramid with the Zr atom on the principal axis. For n = 7-11, their ground state structures can be regarded as attaching one, two, three, four, and five Si atoms to the face of the ZrSi6 -cluster, respectively. That is, the distorted pentagonal bipyramid of ZrSi6 -is the basic structural unit for n = 6-11. For n = 12-16, their global minima are encapsulated shapes with the Zr atom residing in the Si cage. The HP, capped HP, THSQ (three hexagon and six quadrangle), and fullerene cages are evaluated to be the most stable structures of ZrSi12 -, ZrSi13 -, ZrSi14 -, and ZrSi16 -, respectively. For ZrSi15 -, the TPTQ (15m1) and fivecapped PH (pentagonal prism) (15m2) are competitive for the global minimum due to the fact that they are degenerate in energy (the energy difference is only 0.01 eV). Although it is difficult to determine the ground state structure by energy in this situation, the 15m2 isomer is the most probable ground state structure based on simulated PES (see below). To our knowledge, only ZrSi16 -was studied previously by Wu et al [33].
The ground state structures for all di-anionic ZrSin 2-(n = 6-16) are evaluated to be singlet. The growth model of di-anionic ZrSin 2-species is consistent with that of mono-anionic clusters. For n = 6-13, 15, 16, the global minima of di-anionic ZrSin 2-are similar to those of their mono-anions. The cage of di-capped OHFP (one hexagon and four pentagon) is calculated to be the ground state structure for ZrSi14 2-, which differs from that of the corresponding mono-anion.
The distribution of ZrSin 0/-/2-clusters at room-temperature by Boltzmann statistics can further The ground state for all mono-anionic ZrSin -(n = 6-16) is calculated to be a doublet. Their structural growth pattern prefers the attaching type to the encapsulated shape with the Zr atom centered in the Si cage. The global minimum of ZrSi6 -is a distorted pentagonal bipyramid with the Zr atom on the principal axis. For n = 7-11, their ground state structures can be regarded as attaching one, two, three, four, and five Si atoms to the face of the ZrSi6 -cluster, respectively. That is, the distorted pentagonal bipyramid of ZrSi6 -is the basic structural unit for n = 6-11. For n = 12-16, their global minima are encapsulated shapes with the Zr atom residing in the Si cage. The HP, capped HP, THSQ (three hexagon and six quadrangle), and fullerene cages are evaluated to be the most stable structures of ZrSi12 -, ZrSi13 -, ZrSi14 -, and ZrSi16 -, respectively. For ZrSi15 -, the TPTQ (15m1) and fivecapped PH (pentagonal prism) (15m2) are competitive for the global minimum due to the fact that they are degenerate in energy (the energy difference is only 0.01 eV). Although it is difficult to determine the ground state structure by energy in this situation, the 15m2 isomer is the most probable ground state structure based on simulated PES (see below). To our knowledge, only ZrSi16 -was studied previously by Wu et al [33].
The ground state structures for all di-anionic ZrSin 2-(n = 6-16) are evaluated to be singlet. The growth model of di-anionic ZrSin 2-species is consistent with that of mono-anionic clusters. For n = 6-13, 15, 16, the global minima of di-anionic ZrSin 2-are similar to those of their mono-anions. The cage of di-capped OHFP (one hexagon and four pentagon) is calculated to be the ground state structure for ZrSi14 2-, which differs from that of the corresponding mono-anion.
The distribution of ZrSin 0/-/2-clusters at room-temperature by Boltzmann statistics can further support our results. From Table S1, we can see that the lowest-energy structures account for the largest proportion. This verifies the accuracy of the above ground state structures.

PES of Mono-anionic Species
PES is one of the most important experimental tools to get insight into, and to extract electronic fingerprints from, a variety of atomic and molecular clusters as well as condensed-matter systems. Accordingly, the validity of the predicted ground state structures can be tested by way of comparing their theoretical and experimental PES spectra, in which two norms are used. One is the amount of distinct peaks and their position in the PES spectra, and another is the first VDE and/or AEA. The simulated PES of the global minimal structures coupled with the experimental PES spectra are pictured in Figure 4. The theoretical predicted AEAs and the first VDEs are provided in Table 2 along with experimental data. It can be seen from the simulated PES of ZrSi6 -in the range of ≤5.8 eV that there are four distinct peaks (X, A-C) residing at 3.32, 3.63, 4.33, and 5.37 eV, of which the first three peaks are in excellent agreement with experimental data of 3.40, 3.60, and 4.30 eV [23]. The simulated PES of ZrSi7 -has five different peaks (X, A-D) located at 2.86, 3.49, 4.16, 4.85, and 5.59 eV, of which the first three peaks are in excellent accord with the experimental values of 2.80, 3.20, and 4.00 eV [23]. For ZrSi12 -and ZrSi16 -, there are four discrete peaks (X, A-C) centered at 3.99, 4.30, 5.07, and 5.52 eV, and 2.65, 4.11, 4.70, and 5.61 eV, which reproduce well the experimental data of 4.10, 4.40, 4.90, and 5.50 eV, and 2.82, 4.30, 4.80, and 5.50 eV, respectively [22,23]. For ZrSi10 -, four distinct peaks (X, A-C) residing at 3.35, 3.88, 4.45, and 5.58 eV are obtained. In addition to the first peak (X), the remaining three peaks (A-C) are very close to experimental data of 3.80, 4.90, and 5.60 eV, respectively [23]. For ZrSi14 -, four peaks (X, A-C) centered at 3.18, 3.64, 4.59, and 5.16 eV are also obtained, but only later two peaks appear in the experimental PES and agree well with experimental data of 4.50 and 5.20 eV [23]. For n = 8, 9, 11, and 13, their simulated PES has many different peaks. Unfortunately, only three peaks are observed in their experimental PES. And it seems experimentally that their first peaks are not observed analogous to ZrSi10 -and ZrSi14 -. For n = 15, the simulated PES of 15m1 and 15m2 isomers has four peaks (X, A-C) located at 3.00, 3.93, 4.40, and 4.96 eV, and 3.97, 4.51, 4.94, and 5.62 eV, respectively. Apart from first peak (X), the rest of peaks (A-C) for 15m2 are close to the experimental data of 4.70, 5.10, and 5.50 eV, respectively [23]. On the other hand, the three peaks of the experimental spectra of ZrSi15 -are also simulated by the calculations with a shift of the first three peaks (X, A, and B) of the 15m2 isomer by 0.7 eV to the higher binding energy [23]. So we suggest that the 15m2 geometries are the most probable ground state structure. Comparing their theoretical and experimental AEAs (see Table 2), we can see that the AEAs of ZrSin -(n = 6-16) excluding ZrSi11 -, ZrSi12 -, and ZrSi15 -are agreement with those of experimental data. The mean absolute deviations from experimental values are 0.08 eV, and the largest errors are those of ZrSi7 -and ZrSi13 -which are off by 0.11 eV. In short, good agreement between the theoretical and experimental PES spectra sheds further light on the validity of the predicted most stable structures. In light of our reliable theoretical predictions, we suggest that the experimental PES spectra of ZrSi11 -   The ground state of all neutral ZrSi n (n = 6-16) is predicted to be a singlet. The growth pattern adopts the attaching type to the half-cage motif, and finally to the Zr-encapsulated configuration with the Zr atom centered in a Si cage. The global minimum is a face-capped tetragonal bipyramid for ZrSi 6 . For n = 7-9, their ground state structure can be viewed as attaching one, two, and three Si atoms to the face of the ZrSi 6 cluster, respectively. For n = 10-13, they are half-cage motifs consisting of two five-membered rings, a seven-membered and four-membered ring, a seven-membered and five-membered ring, and a seven-membered and six-membered ring, respectively. For n = 14-16, their ground state structures are encapsulated configurations with the Zr atom located in the Si cage. The cage of ZrSi 14 is a distorted HP with two Si atoms symmetrically capping the lateral prism faces. The cage comprised of two pentagons and ten quadrilaterals (TPTQ) is predicted to be the global minimum for ZrSi 15 . A fullerene cage is calculated to be the ground state structure for ZrSi 16 . It is noted that the ground state structures for n = 7-15 presented herein are different from those reported in [18]. In addition, the most stable structures for n = 8-13 differ from those presented in [14], and for n = 14 differs from those presented in [15,17]. The ground state for all mono-anionic ZrSi n -(n = 6-16) is calculated to be a doublet. . Although it is difficult to determine the ground state structure by energy in this situation, the 15m2 isomer is the most probable ground state structure based on simulated PES (see below). To our knowledge, only ZrSi 16 was studied previously by Wu et al. [33]. The ground state structures for all di-anionic ZrSi n 2-(n = 6-16) are evaluated to be singlet.
The growth model of di-anionic ZrSi n 2species is consistent with that of mono-anionic clusters.
The cage of di-capped OHFP (one hexagon and four pentagon) is calculated to be the ground state structure for ZrSi 14 2-, which differs from that of the corresponding mono-anion.
The distribution of ZrSi n 0/-/2clusters at room-temperature by Boltzmann statistics can further support our results. From Table S1, we can see that the lowest-energy structures account for the largest proportion. This verifies the accuracy of the above ground state structures.

PES of Mono-Anionic Species
PES is one of the most important experimental tools to get insight into, and to extract electronic fingerprints from, a variety of atomic and molecular clusters as well as condensed-matter systems. Accordingly, the validity of the predicted ground state structures can be tested by way of comparing their theoretical and experimental PES spectra, in which two norms are used. One is the amount of distinct peaks and their position in the PES spectra, and another is the first VDE and/or AEA. The simulated PES of the global minimal structures coupled with the experimental PES spectra are pictured in Figure 4. The theoretical predicted AEAs and the first VDEs are provided in Table 2  Apart from first peak (X), the rest of peaks (A-C) for 15m2 are close to the experimental data of 4.70, 5.10, and 5.50 eV, respectively [23]. On the other hand, the three peaks of the experimental spectra of ZrSi 15 are also simulated by the calculations with a shift of the first three peaks (X, A, and B) of the 15m2 isomer by 0.7 eV to the higher binding energy [23]. So we suggest that the 15m2 geometries are the most probable ground state structure. Comparing their theoretical and experimental AEAs (see Table 2 In short, good agreement between the theoretical and experimental PES spectra sheds further light on the validity of the predicted most stable structures. In light of our reliable theoretical predictions, we suggest that the experimental PES spectra of ZrSi 11 and ZrSi 15 species should be checked further.
The cause lies in that their theoretical and experimental AEAs and VDEs have large differences. and ZrSi15 -species should be checked further. The cause lies in that their theoretical and experimental AEAs and VDEs have large differences.

Stability
The  The inserts show the experimental PES spectra taken from [22,23].

Stability
The From Figure 5a we can see the data for E b : di-anionic compounds > mono-anionic species > neutral clusters. The higher E b , the more stable the relative stability. It is to say that Zintl anionic Zr-doped Si clusters are more stable than their neutral counterparts, particularly for di-anionic compounds. The cause lies in that neutral ZrSi n compounds still have a dangling bond. When they obtain electrons, these extra electrons decrease the dangling bonds and improve the stability; (Figure 5b) for neutral ZrSi n with n = 9, 12, and 14, for mono-anionic ZrSi n with n = 8, 11, and 13, and for Zintl di-anionic ZrSi n 2with n = 7, 10, 12 and 15, the compounds are more stable than their neighbors. These results are obviously reproduced in Figure 5b because the ∆ 2 E is a sensitive estimate for relative stability. From Figure 5a we can see the data for Eb: di-anionic compounds > mono-anionic species > neutral clusters. The higher Eb, the more stable the relative stability. It is to say that Zintl anionic Zrdoped Si clusters are more stable than their neutral counterparts, particularly for di-anionic compounds. The cause lies in that neutral ZrSin compounds still have a dangling bond. When they obtain electrons, these extra electrons decrease the dangling bonds and improve the stability; ( Figure  5b) for neutral ZrSin with n = 9, 12, and 14, for mono-anionic ZrSinwith n = 8, 11, and 13, and for Zintl di-anionic ZrSin 2-with n = 7, 10, 12 and 15, the compounds are more stable than their neighbors. These results are obviously reproduced in Figure 5b because the ∆ 2 E is a sensitive estimate for relative stability. Egap as a significant physical parameter can be regarded as an index of chemical reactivity, particularly for photochemical reactivity. Larger values of the Egap stand for weaker chemical reactivity. Baerends et al. found that the Egap predicted by pure DFT is closer to the real optical gap than that calculated by hybrid DFT [34]. In light of this finding, the PBE Egap of ZrSin 0/-/2-(n = 6-16) compounds are pictured in Figure 5c as functions of the size of Si clusters, and mPW2PLYP Egap is pictured in Figure S4 in Supplementary Information. It can be seen from Figure S4 that the Egap of mPW2PLYP is larger than that of PBE by 2.46 eV. The Egap of PBE is a very good approximation to the optical gap. The PBE Egap of 1.56 eV for ZrSi16 is in good agreement with an approximate experimental value of 1.36 eV [22]. For n = 6-9, the di-anionic Egap is very small and ranges from 0.27 to 0.81 eV, which is narrower than that of corresponding neutral and mono-anionic compounds. The di-and mono-anionic Egap for n = 10-13, 15, and 16 differ little from each other. For n = 14, the mono-anionic Egap is smaller than that of the di-anion. The neutral Egap for n = 10 and 14 is very close to that of corresponding di-anion. The neutral Egap for n = 11 and 16 is larger than that of di-anion. The neutral Egap for n = 12, 13, and 15 is smaller than that of di-anion. The Egap of 2.25 eV for Zintl di-anionic ZrSi15 2- E gap as a significant physical parameter can be regarded as an index of chemical reactivity, particularly for photochemical reactivity. Larger values of the E gap stand for weaker chemical reactivity. Baerends et al. found that the E gap predicted by pure DFT is closer to the real optical gap than that calculated by hybrid DFT [34]. In light of this finding, the PBE E gap of ZrSi n 0/-/2-(n = 6-16) compounds are pictured in Figure 5c as functions of the size of Si clusters, and mPW2PLYP E gap is pictured in Figure S4 in Supplementary Information. It can be seen from Figure S4 that the E gap of mPW2PLYP is larger than that of PBE by 2.46 eV. The E gap of PBE is a very good approximation to the optical gap.
The PBE E gap of 1.56 eV for ZrSi 16 is in good agreement with an approximate experimental value of 1.36 eV [22]. For n = 6-9, the di-anionic E gap is very small and ranges from 0.27 to 0.81 eV, which is narrower than that of corresponding neutral and mono-anionic compounds. The di-and mono-anionic E gap for n = 10-13, 15, and 16 differ little from each other. For n = 14, the mono-anionic E gap is smaller than that of the di-anion. The neutral E gap for n = 10 and 14 is very close to that of corresponding di-anion. The neutral E gap for n = 11 and 16 is larger than that of di-anion. The neutral E gap for n = 12, 13, and 15 is smaller than that of di-anion. The E gap of 2.25 eV for Zintl di-anionic ZrSi 15 2is the largest among all these compounds. That is, the Zintl di-anionic ZrSi 15 2cluster possesses not only thermodynamic stability, but also chemical stability, which may make it the most suitable building block for novel multi-functional nanomaterials.

Chemical Bonding Analysis
In order to garner an in-depth comprehension of the ideal chemical and thermodynamic stability of the Zintl di-anionic ZrSi 15 2five-capped pentagonal prism, the character of the bonding between the Zr atom and the Si 15 2cage is examined by the AdNDP scheme, which is a generalized NBO (Natural Bond Orbital) search strategy to analyze the delocalized and localized multicenter bonds (encoded as nc-2e, where n can range from 1 (lone pair) to the number of atoms in the cluster) introduced by Zubarev and Boldyrev [35]. As can be seen from Figure 6, the chemical bonding of 66 valence electrons can be split into five types: lone-pair, 2c-2e, 3c-2e, 4c-2e, 10c-2e, and 16c-2e. Each capped Si atom has a lone pair. The pentagonal prism is characterized by ten 2c-2e localized Si-Si σ bonds with 1.86 |e| in each bond. The ten 3c-2e, five 4c-2e, two 10c-2e, and 16c-2e delocalized π+σ bonds with 1.79-2.00 |e| in each bonds are accountable for the interplay between the inner Zr atom and the outer shell of Si 15 and enhance the stability of the encapsulated ZrSi 15 2molecule. is the largest among all these compounds. That is, the Zintl di-anionic ZrSi15 2-cluster possesses not only thermodynamic stability, but also chemical stability, which may make it the most suitable building block for novel multi-functional nanomaterials.

Chemical Bonding Analysis
In order to garner an in-depth comprehension of the ideal chemical and thermodynamic stability of the Zintl di-anionic ZrSi15 2-five-capped pentagonal prism, the character of the bonding between the Zr atom and the Si15 2-cage is examined by the AdNDP scheme, which is a generalized NBO (Natural Bond Orbital) search strategy to analyze the delocalized and localized multicenter bonds (encoded as nc-2e, where n can range from 1 (lone pair) to the number of atoms in the cluster) introduced by Zubarev and Boldyrev [35]. As can be seen from Figure 6, the chemical bonding of 66 valence electrons can be split into five types: lone-pair, 2c-2e, 3c-2e, 4c-2e, 10c-2e, and 16c-2e. Each capped Si atom has a lone pair. The pentagonal prism is characterized by ten 2c-2e localized Si-Si σ bonds with 1.86 |e| in each bond. The ten 3c-2e, five 4c-2e, two 10c-2e, and 16c-2e delocalized π+σ bonds with 1.79-2.00 |e| in each bonds are accountable for the interplay between the inner Zr atom and the outer shell of Si15 and enhance the stability of the encapsulated ZrSi15 2-molecule.

Conclusions
The systematic isomer search for low energy structures of neutral, mono-, and di-anionic Zrdoped Si clusters ZrSin 0/-2-(n = 6-16) are executed by employing an ABCluster global search technique combined with a mPW2PLYP double hybrid density functional scheme. In terms of the evaluated energies, AEA, VDE, and agreement between simulated and experimental PES, the true global minimal structures are confirmed. The results reveal that structural evolution patterns for neutral ZrSin clusters prefer the attaching type (n = 6-9) to the half-cage motif (n = 10-13), and finally to the Zr-encapsulated configuration with the Zr atom centered in a Si cage (n = [14][15][16]. For Zintl mono-and di-anionic ZrSin -/2-, their growth patterns adopt the attaching configuration (n = 6-11) to the encapsulated shape (n = 12-16). The calculated AEA for ZrSin -(n = 6-16) with the exception of n = 11, 12, and 15 are in excellent accord with those of experimental data. The mean absolute deviations from experimental values are 0.08 eV. The further analyses of stability and chemical bonding make it known that two extra electrons not only perfect the structure of ZrSi15 but also improve its chemical and thermodynamic stability, which may make it the most suitable building block for novel multifunctional nanomaterials. We think that this study will provide strong motivation for further

Conclusions
The systematic isomer search for low energy structures of neutral, mono-, and di-anionic Zr-doped Si clusters ZrSi n 0/-2-(n = 6-16) are executed by employing an ABCluster global search technique combined with a mPW2PLYP double hybrid density functional scheme. In terms of the evaluated energies, AEA, VDE, and agreement between simulated and experimental PES, the true global minimal structures are confirmed. The results reveal that structural evolution patterns for neutral ZrSi n clusters prefer the attaching type (n = 6-9) to the half-cage motif (n = 10-13), and finally to the Zr-encapsulated configuration with the Zr atom centered in a Si cage (n = 14-16). For Zintl mono-and di-anionic ZrSi n -/2-, their growth patterns adopt the attaching configuration (n = 6-11) to the encapsulated shape (n = 12-16). The calculated AEA for ZrSi n -(n = 6-16) with the exception of n = 11, 12, and 15 are in excellent accord with those of experimental data. The mean absolute deviations from experimental values are 0.08 eV. The further analyses of stability and chemical bonding make it known that two extra electrons not only perfect the structure of ZrSi 15 but also improve its chemical and thermodynamic stability, which may make it the most suitable building block for novel multi-functional nanomaterials. We think that this study will provide strong motivation for further verification of experimental photoelectron spectroscopy of ZrSi 11 and ZrSi 15 and for further experimental and theoretical studies of other transition metal-doped Si clusters.

Theoretical Methods
We utilized three schemes to search the initial isomers of neutral and Zintl anionic ZrSi n 0/-/2-(n = 6-16) in order to search for their true global minimum. First, the calculations started from unbiased searches for the low-lying structures of series ZrSi n 0/-(n = 6-16) (with the exception of ZrSi n 2-, of which isomers were chosen in light of ZrSi n -) through the ABCluster global search scheme coupled with the GAUSSIAN 09 software suite [36,37]. At least 300 configurations produced by the ABCluster algorithm for each ZrSi n 0/-(n = 6-16) species were optimized at the PBE [38] level with a cc-pVDZ-PP basis set for Zr atoms and a 6-31G basis set for Si atoms [39]. Second, the "substitutional structure", which is derived from the ground state structure of Si n+1 by replacing a Si atom with a Zr atom, was utilized. Third, the structures reported in the preceding publications were utilized [9][10][11][12][13][14][15][16][17][18]. Then, the selected low-lying structural candidates were reoptimized at the PBE level with the MIDDLE (cc-pVTZ basis set for Si atoms and basis set for Zr atoms unchanged) basis sets. Vibrational frequency evaluations were also executed at the same level to make sure all the optimized isomers were true local minimum structures. After accomplishment of the preliminary structure optimization through the PBE scheme, the low-lying structural candidates were, again, selected and reoptimized by using the DH-DFT of mPW2PLYP with the MIDDLE basis set [40]. The mPW2PLYP frequency was not calculated.
To further refine the energies, single-point energy calculations were finally executed at the mPW2PLYP level with the LARGE (aug-cc-pVTZ basis set for Si atoms and basis set for Zr atoms unchanged) basis set. Moreover, NPA were also executed to further comprehend the interplay between the Zr atom and Si clusters. The PES spectra of the negatively charged ions ZrSi n -(6-16) were simulated in terms of Koopman's theorem at the mPW2PLYP/LARGE level via the Multiwfn software suite and compared with the experimental ones [41,42]. The spin multiplicity of singlets, triplets, and quintuplets were considered for neutral ZrSi n and di-anionic ZrSi n 2-(n = 6-16), and doublets and quartets were considered for negatively charged ZrSi n -(n = 6-16) clusters. To check the reliability of the calculations, the ZrSi structures and properties predicted by mPW2PLYP and CCSD(T) were compared and provided in Table 3. From Table 3 we can see that the order of increasing energy predicted by mPW2PLYP is singlet, to quintet, to the triplet state, which is in agreement with the result evaluated by CCSD(T), while B3LYP is not. On the other hand, the bond dissociation energy (D 0 ), AEA, and VDE calculated by the mPW2PLYP method were 2.74, 1.44, and 1.45 eV, respectively, which are in good accordance with experimental values of 2.95 [25], 1.46 [24], and 1.58 eV [24]. In particular, the D 0 and AEA predicted by mPW2PLYP among these methods listed in Table 3 is the closest to the experimental data. As a matter of fact, selection of calculation method and basis set is hardly a trivial matter. The cause lies in that the structures and properties predicted by different methods and basis sets may be different. There are many such example. For instance, at the CCSD(T) level, the LanL2DZ basis set for Si atom predicted the ground state of ZrSi was to be triplet state, but 6-311+G(3df,3pd) basis set for Si atom predicted to be singlet state [1].  Figure S1. Low-lying isomers of neutral ZrSi n (n = 6-16) clusters, point group and relative energy (in eV). Figure S2. Low-lying isomers of mono-anionic ZrSi n -(n = 6-16) clusters, point group and relative energy (in eV). Figure S3. Low-lying isomers of di-anionic ZrSi n 2-(n = 6-16) clusters, point group and relative energy (in eV). Figure S4. The HOMO-LUMO energy gap (E gap ) of ZrSi n 0/-/2-(n = 6-16) clusters. (m) stands for mPW2PLYP calculations and (p) stands for PBEPBE calculations. Table S1. Conformational population (%) for low-lying geometries of ZrSi n 0/-/2 species (n = 6-16) clusters.

Conflicts of Interest:
The authors declare no conflict of interest.