Insights into the structures and electronic properties of Cun+1μ and CunSμ (n = 1–12; μ = 0, ±1) clusters

The stability and reactivity of clusters are closely related to their valence electronic configuration. Doping is a most efficient method to modify the electronic configuration and properties of a cluster. Considering that Cu and S posses one and six valence electrons, respectively, the S doped Cu clusters with even number of valence electrons are expected to be more stable than those with odd number of electrons. By using the swarm intelligence based CALYPSO method on crystal structural prediction, we have explored the structures of neutral and charged Cun+1 and CunS (n = 1–12) clusters. The electronic properties of the lowest energy structures have been investigated systemically by first-principles calculations with density functional theory. The results showed that the clusters with a valence count of 2, 8 and 12 appear to be magic numbers with enhanced stability. In addition, several geometry-related-properties have been discussed and compared with those results available in the literature.

As we know that the structures and properties of clusters are very sensitive to the number of atoms, which can change dramatically with the addition or substitution of one atom. The introduction of a doped atom in copper clusters, such as a sulfur atom in our work, can undoubtedly change the clusters' structure, which in further alters their physical and chemical properties significantly. Copper sulfides have attracted sustained research interests over the past decades due to their wide applications in solar cell devices, nonlinear optical materials, lithium ion batteries, nanometer-scale switches and gas sensors [22][23][24][25] . For example, the electronic structures of Cu 2 S cluster were calculated in the local density functional approximation 26 . Their theoretical results suggested that the structures closely resemble an S − bridging a Cu 2 + as shown by the geometries and charge distributions. The potential energy surface minima of neutral copper sulfide clusters have been explored using the Coalescence Kick global optimization method and combined with a DFT approach 27 . The calculated HOMO-LUMO gap (1.3 to 3.3 eV) showed that the (CuS) n clusters can be considered as a suitable candidate for renewable energy sources. At last, Scott et al. 's results confirmed that Cu-S-Cu units can be found in an active site of some metalloproteinase like cytochrome c-oxidase 28 . The short Cu-Cu distance and small Cu-S-Cu bond angle play an essential role in the electron transport made by this protein.
It is well known that the Jellium model plays a significant role in the realm of cluster. It incorporates a positively charged (spherical or non-spherical) background potential, which results in discrete energy levels of the delocalized ("metallic") electrons corresponding to angular momentum shells (in the spherical case the shells are labeled as: 1S 2 1P 6 1D 10 2S 2 1F 14 , 2P 6 , 1G 18 …) 29 . There are certain valence electronic configurations (2,8,18,20,40, 58…) know as magic numbers, that exhibit increased stability relative to their neighboring configurations. So, these clusters are generally resistant to reactivity with small clusters [30][31][32][33][34] . For example, Jaque et al. 14 studied the molecular structures, binding energy, electronic properties and reactivity descriptors for the neutral copper cluster, and found that the large HOMO-LUMO energies gap, minimum polarizability and maximum hardness corresponding with the 2 and 8 valence electrons are operative for characterizing and rationalizing the electronic properties of copper clusters. However, other groups found that some clusters with the non-magic number total electrons can also exhibit an increased stability 35,36 . For example, Rebere et al. 36 reported that Ag 15 + , Ag 14 and Ag 13 − clusters with 14 valence electrons are resistant to reactivity with oxygen, even though they do not have a magic number of electrons. This simulated further and more detailed investigations on the magic numbers.
Although the structural and electronic properties of copper and sulfur doped copper clusters have been investigated by both experiment and theory, some important question arises: (i) For different copper clusters with the small energy difference between the atomic s and d levels, how does hybridization effect from different molecular orbitals? (ii) How does the structural vary as a consequence in anionic and cationic Cu n S clusters? Furthermore, are their structures and properties greatly distinct from the neutral Cu n S clusters? (iii) In view of the fact that the electronic structure of Cu and Ag featuring a closed d shell and a singly occupied s shell (3d 10 4s 1 and 4d 10 5s 1 ), whether the Jellium model is in accordance with the pure copper clusters? With a doped sulfur, how does the magic number vary? (iv) For the charged Cu n S clusters, how does HOMO-LUMO gap change? Hence it is necessary to carry out a systematic investigation to reveal the magic number, structures and electronic properties of copper and sulfur doped copper clusters.
With this purpose in mind, we perform a systematic, in-depth study of the geometric structures and electronic properties for neutral Cu n+1 and Cu n S (n = 1-12) clusters and their ions based on particle swarm optimization algorithm and density function theory. Thus, we are able to advance a fundamental understanding of the ground state geometric configurations. In addition, the magic numbers, HOMO-LUMO gaps, density of states (DOS), adaptive natural density partitioning (AdNDP), electron localization function (ELF) and Mayer bond order are also analyzed. Subsequently, the calculated values are compared with available experimental and theoretical data.

Results and Discussions
Structural properties. We performed unbiased searches for the global minima structures of neutral and charged Cu n+1 and Cu n S (n = 1-12) clusters. Previously reported structures are successfully reproduced with CALYPSO method from experimentally and theoretically. The global minimum structures for each size are displayed in Fig. 1. The corresponding electronic states, symmetries, total energies are also determined and listed in Table S1 in the Supporting Information (S1).
As depicted in Fig. 1, it is clear that the charges have important influence on the structures of pure and doped clusters. We found nearly the same structure only in the cases of n = 1, 3 for Cu n+1 and n = 1, 2 for Cu n S, respectively, while the remaining clusters contain different structures. For example, both Cu 3 + and Cu 3 have a triangular structure, while Cu 3 − is composed of a linear structure. The same trend can be also observed for Cu 3 S +/0/− cluster, namely, Cu 3 S + and Cu 3 S possess the same triangular structure, whereas Cu 3 S − has a rhombus structure. In addition, the transition from two-dimensional (2D) to three-dimensional (3D) configurations is closely related to their charge state. For non doped-cluster, the 2D-3D structure transition occurs at n = 4-5 for Cu n + , n = 6-7 for Cu n , and n = 5-6 for Cu n − ; while for S doped ones, n = 2-3 for Cu n S + , n = 2-3 for Cu n S, and n = 5-6 for Cu n S − , respectively. For neutral copper cluster, two growth patterns can be obtained here. One is to remain essentially two dimensional and grow by forming successive triangle faces for n = 2-6. From n > 7, the cluster forms a pentagon bipyramid structure with the two capping atoms forming a bond. And, it becomes more closely packed with the increasing atomic number of clusters. For charged copper cluster, the ground state structures exhibit a layer-like 3D configuration from n = 5 to 6, which display a preference for layered and pyramidal geometries. In addition, our structures for Cu 2  20 . Nevertheless, it should be mentioned that a new lowest energy structure of Cu 9 cluster is identified, which is different from previous studies 14- 16 . In present work, a C 2v symmetry Cu 9 cluster is found as the lowest energy structure; however, different isomer with C s symmetry are reported by Jug (ALLCHEM), Jaque (B3PW91/LANL2DZ) and Ramirez et al.
(BLYP/6-311+G(d)). The different lowest energy structures may be due to the fact that the different functional and basis sets may reverse the order of clusters in some cases. So, in order to confirm the lowest energy structure of Cu 9 cluster, the same functional and basis sets (B3PW91/LANL2DZ) are performed based on the DFT calculations. Results showed that the Cu 9 structure C s , reported in ref. 14 as a minimum, in our calculation is found to be the second one, at a relative energy of 0.02 eV above our minimum. In addition, it is also noted that the Gibbs free energy of Jaque et al. 's results is 6.3 × 10 3 J/mol higher than our calculated work. Above analysis indicated that that our results are energetically lower and structurally more stable than Jaque et al. 's results.
For the neutral and charged Cu n S clusters, different types of clusters possess the same growth pattern with rare exceptions. Practically speaking, the doped clusters are formed by adding a copper atom into the smaller sized Cu n−1 S clusters. Any additional S atom is just one more vertex to surround the Cu atom, and finally the cluster forms a triangular prism. In conclusion, the growth pattern of Cu n S +/0/− clusters are quite different from Cu n +/0/− clusters. One Cu atom directly added on the Cu n−1 S clusters are dominant growth pattern for Cu n S +/0/− clusters.
To check the accuracy of the lowest energy structures, several available experimental and theoretical results are calculated and compared. For Cu n clusters, the evolution of VIP, VEA and static mean polarizablity (α) values is plotted in Figure S1, together with available experimental and calculated results 3-5, 13, 16 . For Cu n − clusters, the calculated ADE and VDE are listed in Table 1, which also includes the available experimental values taken from ref. 1. For Cu 2 S cluster, the calculated bond length (Cu-S, Cu-Cu) and bond angle (CuSCu) (2.14 Å, 2.67 Å and 77.5°) are close to the Boris and Mahe et. al's. work (2.11 Å, 2.63 Å and 76.7°) 37,38 . In addition, we found that our  calculated results about the bond length and angles, dissociation energies and frequency are also in excellent agree with the experimental and theoretical results for Cu n S clusters 26,27,37,38 . In conclusion, we can see that our calculated values are in good agreement with the experimental and theoretical results. Thus, we believe that the choice of functional and basis sets could be reasonably good to describe the present systems.
Stabilities and electronic properties. To analyze the relative stabilities of the lowest energy structures of Cu n+1 +/0/− and Cu n S +/0/− (n = 1-12) clusters, the averaged binding energies E b (n) and second-order difference of energies ∆ 2 E(n) have been derived. E b (n) and ∆ 2 E(n) of each complex for Cu n S +/0/− clusters are calculated as follows: For Cu n +/0/− clusters, E b (n + 1) and ∆ 2 E(n + 1) are defined as: are relatively more stable than its neighboring clusters. For Cu n S −/+ clusters, Cu 3 S − , Cu 5 S − , Cu 7 S − , Cu 9 S − and Cu 3 S + have a higher relative stability than their neighbors, respectively. (iii) For Cu n+1 +/0/− clusters, E b (n) have a monotonically increasing with an increase of copper atoms, revealing an enhanced effect on the stabilities of the Cu n+1 +/0/− clusters with the cluster size increasing. In addition, the charged clusters have larger values for the binding energy than neutral ones, while ionic clusters possessing the largest values lie in between them. Another important parameter to evaluate the relative stabilities is the size-dependent second energy difference, the values of neutral and charged clusters exhibit obvious odd-even alternations as shown in Fig. 2. In which the charged and neutral clusters are more stable with even and odd cluster size n for non-doped Cu n+1 +/0/− , respectively. Whereas the charged and neutral clusters are more stable with odd and even cluster size n for doped Cu n S +/0/− , respectively. The above analysis can be explained by this can be attributed to the presence or absence of unpaired electrons of the complex. In addition, we can easily point to some conspicuous peaks: n = 2, 7, 6 for Cu n+1 +/0/− clusters as well as n = 3, 2, 5 for Cu n S +/0/− clusters. It suggests that these clusters are more stable compared with other clusters.
The HOMO-LUMO energy gaps (E gap ) have been proved to be a powerful tool to represent the ability of the molecule to participate in the chemical reaction in some degree. Larger values of values indicate stronger chemical stability. The calculated of HOMO-LUMO energy gaps are shown in Fig. 3 for the lowest energy structures of Cu n+1 +/0/− and Cu n S +/0/− (n = 1-12) clusters. It is interesting to notice that Cu n+1 with even number of electrons are more stable than other clusters with odd number of electrons. Contrary, for Cu n+1 + and Cu n+1 − clusters, odd-even alteration behaviors of E gap curve present the opposite trend compared to the corresponding neutral clusters. This can be explained that even numbers electrons are always exist in pairs and can lead to a closed shell electronic structure. We also observe that E gap of cationic copper clusters are larger than their anionic counterparts, indicating that Cu n+1 + clusters are less stable, which is in excellent agreement with the results of averaged binding energies shown in the Fig. 2. As for Cu n+1 − clusters, four local maxima of E gap are found at n = 2, 6, 8 and 12, respectively, suggesting that Cu 3,7,9,13 − clusters are more stable than their neighbors. The local maxima of the E gap for Cu n+1 and Cu n+1 + clusters indicate that neutral Cu 2, 6, 8 and cationic Cu 3, 5, 9 + clusters possess enhanced relative stability. Similarly, for Cu n S +/− clusters, except Cu 2 S + , the results of E gap present the same odd-even oscillation. Moreover, the E gap of Cu n S − are always lower than their cationic counterparts with the exception of CuS + . Those results indicate that cationic clusters are more stable than corresponding anionic clusters. In the case of Cu n S +/0/− clusters, some local maxima E gap values are observed at Cu 2 S, Cu 4 S, Cu 6 S for neutral clusters, Cu 3 S + , Cu 5 S + , Cu 9 S + for cationic clusters as well as Cu 3 S − , Cu 5 S − , Cu 11 S − for anionic clusters. Combining the above analysis on E b (n), ∆ 2 E(n) and E gap values of Cu n+1 +/0/− and Cu n S +/0/− clusters, we can infer that Cu 8 , Cu 3 + , Cu 7 − , Cu 2 S, Cu 3 S + and Cu 5 S − clusters exert a remarkable chemical stability. At the same time, we checked the valence electrons of above stable clusters, the number is 2 and 8 for Cu 3 + and Cu 8 , Cu 7 − , Cu 2 S, Cu 3 S + respectively, which are magic numbers according to Jellium model. However, our theoretical predictions point out that the clusters of Cu 5 S − is particularly stable despite the cluster has 12 valence electrons and do not conform to the predictions of the Jellium model. Finally, for Cu n S +/0/− clusters, the calculated HOMO-LUMO gaps range from 1.01-4.29 eV, 0.8-2.96 eV, 1.26-2.53 eV, respectively. These are desirable band gaps of semiconductor nanomaterials, from an application point of view, it suitable for application in photocatalysis field (especially for Cu n S + clusters). We expect that studies of the neutral and charged sulfur-copper clusters could provide favorable information in the searching functional design and application of the renewable energy sources.
To understand the nature of the chemical bonding and the formation mechanism of these various Cu-S compounds, the total density of states (TDOS) and partial density of states (PDOS) of Cu 2 S and CuS cluster are shown in Fig. 4. In addition, the Figures of Cu 3 S + and Cu 5 S − clusters are discussed and plotted in Figures S2 and S3 (see S1). From Fig. 4, the band gaps of Cu 2 S are larger than that of CuS, and the TDOS of Fermi level for Cu 2 S clusters is higher than CuS clusters, indicating the metallicity of CuS cluster is weakened due to a Cu ejection. We can also note that the TDOS of Cu 2 S clusters mainly comes from copper atoms in the region from −0.4 to −0.2 a.u. and sulfur atom in the region from 0.05 to 0.15 a.u. As evidenced by the diagram, the HOMO level is dominated primarily by the S-p orbital and partially by the Cu-p, d orbitals, in which Cu-p energy level is lower than Cu-d and S-p orbital energy level is much higher than Cu-d. The contribution from Cu-s and S-s orbital is almost zero. At LUMO case, the level is composed of Cu-s, p, d orbitals and S-s, p orbitals. The major contributions come from Cu-s orbital. Moreover, the orbital energy levels of Cu-p and S-p are higher than that of S-s; the contribution from Cu-d orbital is very small. These results provide a certification for our discussion based on molecular orbital composition analysis.
In order to decipher the chemical bonding of cluster, we performed a detailed analysis based on the localized orbitals resulting from the AdNDP method. AdNDP method, which is developed by Zubarev and Boldyrev 39 , is a very efficient and visual approach to the interpretation of the molecular orbital because it is an extension of the natural bond orbital analysis. It represents the molecular electronic structure in terms of n-center two-electron (nc-2e) bonds, the familiar lone pairs (1c-2e) and localized 2c-2e bonds or delocalized nc-2e bonds (3 ≤ n ≤ total number of atoms in the system). Moreover, the occupation numbers (ON) are expected to be close to the maximum values (ON = 2.00|e|).
When performing AdNDP approach to characterize the chemical bonding in Cu 2 S, the full-filled 3d orbital of Cu are not shown in this figure with ON ranging from 1.96 to 1.99 |e|. Excluding Cu (3d) orbital, Cu 2 S has 8 valence electrons with each copper atom contributing one valence electron and each sulfur atom contributing six valence electrons. The results of Cu 2 S cluster are depicted in Fig. 5. AdNDP analyses find four electrons are localized along the two localized Cu-S two-center two-electron (2c-2e) σ-bonds (ON = 1.972|e|). The remaining 4 electrons contain two lone pairs (LPs) with one s-type (ON = 1.956 |e|) and one p-type (ON = 1.936 |e|) on the sulfur atom. More specifically, per 2c-2e bond possesses the atomic contribution of 39.28%Cu + 60.72%S. In Cu-S 2c-2e bond, Cu-s and p offer 34.22% and 2.67% contribution to the Cu-based orbital, the S-based orbital possess the contributions of S-p (52.57%) and S-s (7.76%). Evidently, the Cu-S 2c-2e σ-bond is mainly dominated by the components from Cu-s and S-p type orbitals.
For Cu 3 S + cluster, the 8 valence pairs are localized as one LPs with one S-s (ON = 1.968 |e|) and three localized Cu-S 2c-2e σ-bonds (ON = 1.981 |e|). The atomic contributions come from 32.61%Cu + 67.39%S for per 2c-2e bond. In Cu-S 2c-2e bond, the Cu-based orbital consist of Cu-s (27.10%) and Cu-p (3.96%); the corresponding S-based orbital is composed of 51.18% S-p and 16.22% S-s. Obviously, Cu-s and S-p possess three localized Cu-S 2c-2e bond. In addition, there is no covalent bonding between Cu and Cu atoms, indicating that interaction between Cu and Cu atoms is a kind of non-Lewis interaction.
For Cu 3 + cluster, the results found that a pair of electrons is located inside of the copper triangle being delocalized over all three copper atoms. The occupation number of 3c-2e Cu-Cu-Cu σ-bonds is equal to the ideal value of 2.00 |e|. Due to the same contribution from the different three copper atoms, the corresponding atomic contribution is 33.33%Cu(1) + 33.33%Cu(2) + 33.33%Cu (3). Moreover, the 3c-2e Cu-Cu-Cu σ-bond possesses the orbital contribution of Cu-s (30.87%) and p (2.46%). The major contribution to the 3c-2e σ-bonds come from three Cu-s type orbital.
In order to study the bond interactions of doped systems, we conducted the electron localization function (ELF) analysis using the Multiwfn program package 40 . The ELF is usually interpreted as the probability of finding an electron pair localized in some region of the real space. There are three special cases: ELF = 0 indicates that there is no electron density between atomic orbitals; ELF = 0.5 shows the regions with bonding of a metallic character; ELF = 1 corresponds to perfect localization, and means covalent bonds and inner shell or lone pair electrons. To elucidate the type of bonding in clusters, bond length and Mayer bond order and ELF cut planes for Cu 2 S, Cu 3 S + and Cu 5 S − cluster are illustrated in Fig. 6. We can see that the ELF values near Cu atoms are very high which associated with 1s core electrons, whereas the ELF values of regions between copper atoms are very small (nearly zero), it means that the chemical bonding between Cu and Cu atoms can be neglected. The ELF values of bonds are small between Cu and S atoms, indicating that there exist typical of ionic bonding in Cu 2 S, Cu 3 S + and Cu 5 S − clusters. In addition, the distribution of ELF value near the S atom is not spherically symmetric, and there is a region with higher localized degree. This implies that the S atom through sp hybridization from normal two-center σ bonds with the near copper atoms. The present results are in good agreement with AdNDP chemical bonding analysis. Moreover, in Cu 5 S − cluster, it is clear that electrons are mainly localized within the two outside triangles, corresponding to the three-center two-electron bond (3c-2e). Simultaneously, Mayer three-center bond order analyses show that the value of outside triangle is 0.071e. In addition, based on the two-center bond order analysis, the bond order of boundary Cu-Cu bond (0.521) is markedly larger than the inner one (0.290). Results indicate that the bond of boundary Cu-Cu is more stable than the inner ones; the bond of boundary Cu-S is more stable than that of Cu-Cu. At last, comparing the different bond lengths, we can find that short bond lengths lead to strong interatomic interactions, which play an important role in stabilizing structures.

Conclusions
In this paper, we have performed a global minimum search for the lowest energy structures of neutral, cationic and anionic Cu n+1 and Cu n S (n = 1-12) clusters by using CALYPSO method in combination with density functional theory. The results are summarized as below: (i) The structure searches show that the lowest energy structures are sensitive to the charge states. A new global minimum state of Cu 9 cluster is obtained, which is more stable than previous work. For Cu n S +/0/− clusters, Cu directly added on the Cu n−1 S clusters to form the Cu n S clusters is the dominant growth pattern. (ii) Sulfur atom doping has considerable influence on not only geometries and properties but also the valence electron count of copper clusters. Trends of the atomic binding energies, second-order difference of energies and HOMO-LUMO gaps showed that the clusters containing even number of electrons maintain greater stability than odd number of electrons. In addition, Cu 3 + , Cu 8 , Cu 7 − , Cu 2 S and Cu 3 S + follow the trends predicted by the Jellium model with the 2, 8 valence electron systems being the most stable. However, Cu 5 S − cluster with 12 valence electrons does not correspond to the magic numbers also exhibit an increased stability. (iii) The calculated HOMO-LUMO gaps range from 1.01-4.29 eV, 0.8-2.96 eV, 1.26-2.53 eV, respectively, which make Cu n S +/0/− clusters suitable candidates for renewable energy sources (especially Cu n S + clusters). Density of states reveals that the orbital contribution of HOMO level is dominated by the S-p and Cu-s, p, d, the orbital contribution from S-s is almost zero. In LUMO level, the orbital contributions are composed of Cus, p, d and S-s, p, respectively. Moreover, the structural symmetry effect corresponding orbital composition. (iv) AdNDP analyses show the chemical bonding patterns are in good agreement with the geometric structures and the ELF and Mayer Bond order analysis. The bond of boundary Cu-Cu is more stable than the inner ones. The bond of boundary Cu-S is more stable than that of Cu-Cu.

Methods
The lowest energy structures of neutral and charged Cu n+1 and Cu n S clusters are searched by the swarm-intelligence based CALYPSO structure prediction method [41][42][43] . This method is based on globally minimizing potential energy surfaces, merging ab initio total energy calculations with CALYPSO cluster prediction through particle swarm optimization. It has been successful in correctly predicting structures for various systems [44][45][46] . In this prediction method, each generation contain 20 structures, 70% of which are generated by particle swarm optimization (PSO). PSO is a population approach based stochastic optimization technique developed by Eberhart and Kennedy in 1995 47,48 . The others are new and will be generated randomly. In process of searching, a sequence of 50 generations of structural candidates is followed to achieve convergence. So, we can achieve 1000 structurally different low-lying isomers. Among the 1000 isomers, the top fifty low-lying isomers are collected as candidates. Those structures with energy difference from the lowest energy isomers less than 0.3 eV are further optimized to identify the lowest energy structure. The further geometry optimizations are performed with no symmetry constraints at the level of the generalized gradient approximation (GGA) using the exchange-correlation functional B3P86 49,50 as implemented in the GAUSSIAN09 package 51 . The basis set labeled GENECP, i. e. the LanL2DZ basis set for Cu atom and 6-311 + G* basis set for S atom, respectively 52,53 .