Simultaneous prediction of atomic structure and stability of nanoclusters in a wide area of compositions

TL;DR

This paper introduces a universal, efficient method for predicting atomic structures and stability of nanoclusters across a wide compositional space, enabling comprehensive first-principles studies.

Contribution

The authors develop a joint evolutionary algorithm that accelerates cluster structure prediction by leveraging structural similarities, significantly improving speed and coverage.

Findings

Speedup of up to 50 times over existing methods

First-principles global optimization of 315 Si-O clusters

Discovery of both expected and novel stable clusters

Abstract

We present a universal method for the large-scale prediction of the atomic structure of clusters. Our algorithm performs the joint evolutionary search for all clusters in a given area of the compositional space and takes advantage of structural similarities frequently observed in clusters of close compositions. The resulting speedup is up to 50 times compared to current methods. This enables the first-principles studies of multi-component clusters with full coverage of a wide range of compositions. As an example, we report an unprecedented first-principles global optimization of 315 SinOm clusters with n<=15 and m<=20. The obtained map of Si-O cluster stability shows the existence of both expected (SiO2)n and unexpected (e.g. Si4O18) stable ("magic") clusters, which can be important for miscellaneous applications.

Figures

Equations

$${\Delta }_{ss}E(\mathbf{n})=E({\mathbf{n}}_s^{+})+E({\mathbf{n}}_s^{-})-2E(\mathbf{n})$$

$$O_i=\frac{{\sum }_j\exp (-(r_{ij}-R_i-R_j)/d)}{{\max }_j\exp (-(r_{ij}-R_i-R_j)/d)}$$

$$E_{\mathrm{diss}}(n,m,k,l)=E(k,l)+E(n-k,m-l)-E(n,m)$$

Full text

\AbstractOn

Skoltech] Skolkovo Institute of Science and Technology, Skolkovo Innovation Center, Nobel St. 3, Moscow 143026, Russia \alsoaffiliation[LPI] P.N. Lebedev Physical Institute of the Russian Academy of Sciences - 119991 Leninskii prosp. 53, Moscow, Russia

Skoltech] Skolkovo Institute of Science and Technology, Skolkovo Innovation Center, Nobel St. 3, Moscow 143026, Russia \alsoaffiliation[LPI] P.N. Lebedev Physical Institute of the Russian Academy of Sciences - 119991 Leninskii prosp. 53, Moscow, Russia

LPI] P.N. Lebedev Physical Institute of the Russian Academy of Sciences - 119991 Leninskii prosp. 53, Moscow, Russia

Skoltech] Skolkovo Institute of Science and Technology, Skolkovo Innovation Center, Nobel St. 3, Moscow 143026, Russia \alsoaffiliation[MIPT] Moscow Institute of Physics and Technology, Dolgoprudny, Moscow Region 141700, Russia \alsoaffiliation[NPU] Northwestern Polytechnical University, Xi’an, Shaanxi 710072, People’s Republic of China

{tocentry}

The unique properties of nanoparticles are extensively used in optoelectronics, photovoltaics, photocatalysis, biomedicine, etc. These properties are closely linked to the atomic structure of particles, what is more explicit in the small particles and nanoclusters 1, 2. Despite the importance of knowing the structure, its experimental determination remains very difficult 3. For this reason the main body of structural information on clusters is obtained via first-principles calculations 4 which were mostly done either for monoatomic clusters or for binary clusters of stoichiometric composition corresponding to the bulk compounds, while clusters of general composition were studied only in few publications5, 6. Such an accent in ab initio research ignores the fact that the chemistry of clusters is much richer than that of solids because of a large share of surface atoms. Multi-component clusters often have stable compositions, which are far from chemical compounds presented in the bulk $x$–$T$ phase diagram. This is of interest not only for basic chemistry of clusters. It significantly increases the scope of candidate nanomaterials for practical applications such as: the development of efficient and affordable catalysts 7, 8 and magnets 9, the investigation of complex processes of nucleation and particle growth 10, 11, 12, etc.

The bottle-neck of first-principles activity in cluster studies is the computational cost of atomic structure determination, which is a global optimization of the total energy among all possible atomic configurations. There are several methods of structure prediction (basin and minima hopping 13, 14, simulated annealing 15, evolutionary algorithm 16, etc.), however they all involve thousands of local optimizations (relaxations) even for finding a structure of one cluster. In the applications mentioned above, the computation of atomic structure and the screening of stability and properties are required in wide regions including hundreds of cluster compositions, therefore such first-principles investigations turn out extremely exhausting. To reduce the computational cost, the global optimization is frequently performed in combination with semiempirical methods of force fields 17, 18. The success greatly depends on the model potential, which is often difficult to make sufficiently accurate. Here we suggest a different approach to this problem, which does not invoke semiempirical potentials at all. Our method predicts all clusters in the whole given area of compositions simultaneously in a highly efficient manner that incorporates exchange of structural information, i.e. learning between clusters of different compositions. The effectiveness of our approach is based on the frequent similarity of structural motifs in clusters of close compositions. We will refer our technique to as variable-composition cluster search in contrary to the previous, fixed-composition approaches.

Our method is derived from the evolutionary algorithm implemented in the USPEX code16, 19, 20, proved to be successful for predicting novel materials21, 22. Briefly, the algorithm is based on the analogy with natural selection: first population of structures (the 1-st generation) is created randomly and locally optimized. Certain percentage of best structures serves as parents for the next generation, which is produced using the so-called variation operators: cross-over, mutation, etc. Then the next population is locally optimized and so on, until the convergence is achieved.

Our technique performs a joint evolutionary search for the whole given area of compositions at once. It required introducing two major innovations. First, the selection of the fittest configurations (selecting best parents) is performed on equal footing for clusters of all compositions. Second, we developed new variation operators to provide the transfer of structural information between clusters of different compositions.

Let’s consider the new selection procedure in more detail. It is grounded in the notion of “magic” clusters. We classify a cluster as magic if a pair of such clusters is stable against the transfer of one atom between them. Giving it a formal description, we denote a composition of cluster with the formula X${}_{n_{1}}$Y${}_{n_{2}}$…Z${}_{n_{k}}$ by a vector $\mathbf{n}=(n_{1},\,n_{2},\ldots,\,n_{k})$. If an atom of sort $s$ is added or removed, the composition becomes $\mathbf{n}_{s}^{\pm}=(n_{1},\,\ldots,\,n_{s}\pm 1,\ldots,\,n_{k})$. The cluster is magic if the second-order differences of the energy:

[TABLE]

are positive for all sorts $s$. Having the (non-regular) set of magic compositions, we build the reference energy surface $E_{\mathrm{ref}}(\mathbf{n})$ as a piecewise linear interpolation of ’magic’ energies over all given area of compositions. A configuration is classified as ’best’ and selected to participate as a parent for the next generation, if its energy falls in the interval $E_{\mathrm{conf}}(\mathbf{n})-E_{\mathrm{ref}}(\mathbf{n})\leq\Delta E_{\mathrm{sel}}$. The interval $\Delta E_{\mathrm{sel}}$ is defined so that the share of the best configurations is $N^{\mathrm{best}}/N\sim 0.6-0.8$. We note that the global thermodynamic stability of clusters is not required, as any cluster system is unstable with respect to growth or coalescence: the only truly thermodynamically stable cluster is the infinite one, i.e. a crystal. This situation is totally different from the variable-stoichiometry prediction of crystalline structures 23, which provides only the thermodynamically stable phases.

Our method uses all variation operators of the standard fixed-composition approach20: (1) creating structures with random point symmetry; (2) permutations of chemically different atoms; (3) softmutation (displacement of atoms along eigenvectors of the softest vibrational modes) and (4) heredity (creating child structure from fragments of two parents). Among these operators only heredity is suitable for structural exchange between clusters of different compositions and was modified accordingly. To further enhance such exchange, three new variation operators are introduced: (5) transmutation (change of chemical identity of randomly selected atoms), (6) removal of one atom from the cluster, and (7) addition of one atom to the cluster.

For operators (6) and (7) a location, where one atom should be removed or added, is of importance. The choice of an atom $i$ to be removed is defined by its effective coordination number $O_{i}$:

[TABLE]

Here $r_{ij}$ is the distance between atoms $i$ and $j$, $R_{i}$ and $R_{j}$ are the covalent radii of atoms $i$ and $j$ depending only on their chemical identities, and $d=0.23\,\text{\AA}$ is the empirically determined parameter24. The $i$-th atom is removed with a probability $p_{i}$ proportional to $\max_{i\in s}[O_{i}]-O_{i}$, where $s$ is the identity of an atom $i$. Thus, the removal of weakly-bound atoms is preferable. The same applies to the choice of an atom to which an additional atom should be attached. Such atom addition enhances coupling between a weakly bound atom and the remainder of cluster that gives the maximum gain in binding energy.

We note that new add/remove atom operators usually do not create new structural motifs, but spread the good ones between different compositions, thus providing the thorough exploration of the low-energy areas of the landscape. Other operators, considering their greater stochasticity, are responsible for a sufficient level of structural diversity. Such an approach provides a balance between scattering of trials for the effective sampling of search space and focusing in the most promising regions which is crucial for any heuristic optimization scheme.

The seven variation operators show varying efficiency for different systems and even at different stages of evolutionary search. To take this into account we developed a control procedure which changes on-the-fly the number of offspring produced by each operator considering their efficiency at previous stages of the search 25. To increase structure diversity, our method uses yet another, very powerful tool, namely, “antiseeds”. It gives energy penalty for structures which are best for too many generations20. This expedient allows one to reliably determine structures of all non-magic clusters belonging to the area of compositions, as well as low-energy isomer structures.

One may recollect that many researchers tried to exploit structural similarity of clusters of close compositions “by hand”10, 26, 27, 28. However, the number of possible configurations is still extremely large. For this reason such manual sampling requires too much effort and often fails. In contrast, our method makes similar sampling fully automatically, integrating it into the general evolutionary process. Combined functioning of all variation operators, new selection scheme, “antiseed” technique and other features give synergistic effect that results in the high performance of our method.

We test our technique on two model systems: Lennard-Jones clusters with 30-60 atoms and $\text{Si}_{n}\text{O}_{m}$ clusters ($n=6-8$ and $m=10-16$) within the semiempirical MNDO approach as implemented in the MOPAC package29. We compare the convergence speed of the new approach with the standard fixed-composition technique implemented in the USPEX code20. Fig. 1 shows the energy deviations from the ground state averaged over all compositions as a function of the number of relaxations for both methods. This integral characteristic of the convergence rate shows a great speedup of the new method ($\sim 5$ times for LJ clusters and up to 50 times for $\text{Si}_{n}\text{O}_{m}$ clusters). A more thorough analysis presented in Supplementary Information (SI) shows high efficiency of new variation operators which produce most of low-energy offspring of the new method.

As the first real application, we chose $\text{Si}_{n}\text{O}_{m}$ clusters owing to their practical importance and great structural diversity (crystalline silica alone has 14 structural forms). The search was performed at the ab initio level within the unprecedentedly large range of compositions ($1\leq n\leq 15$ and $0\leq m\leq 20$ i.e. 315 cluster compositions). We note that the earlier investigations of $\text{Si}_{n}\text{O}_{m}$ were either done for relatively small clusters ($n\leq 7$) 6, 26, 28, 30, 31 or focused on stoichiometric compounds $(\text{SiO})_{n}$ or $(\text{SiO}_{2})_{n}$4, 10, 27, 32, 33. Even for these compounds it is often seen that newer papers report lower-energy structures than the older ones. As an example, we mention recent study of $(\text{SiO})_{n}$ clusters which were constructed “by hand”, joining fragments of $(\text{SiO}_{2})_{n}$ and $\text{Si}_{n}$ clusters together 10. Surprisingly, these structures turned out to be better than those reported in all previous studies using global optimization techniques.

We perform our global optimization combined with density functional calculations within the PAW-PBE approximation implemented in the VASP code 34, 35. The energies of 10 best structures for each composition were refined using the GAUSSIAN code 36 with the B3LYP/6-311+G(2d,p) approach37. Comparison of our results with earlier publications4, 6, 10, 26, 27, 28, 30, 31, 32, 33 showed that 101 optimal structures of $\text{Si}_{n}\text{O}_{m}$ clusters were correctly reproduced by us, 17 better ones were found, and 197 clusters were studied for the first time to our best knowledge (see table 1 in SI). Fig. 2 shows $\text{Si}_{n}\text{O}_{m}$ stable (magic) clusters which are divided into four groups discussed below. The optimal structures of all 315 clusters are given in SI (table 2).

Scanning over wide composition areas reveals trends in cluster structure and related properties. To illustrate this point, we explore the stability patterns of $\text{Si}_{n}\text{O}_{m}$ clusters using two criteria. The first one characterizes the resistance toward the transfer of Si or O atoms between clusters. It calculates second-order differences over Si and O atoms ($\Delta E_{nn}(n,m)$ and $\Delta E_{mm}(n,m)$ of Eq. (1)) and takes the minimal one: $\Delta_{\mathrm{min}}(n,m)=\min\{\Delta_{nn}E(n,m),\Delta_{mm}E(n,m)\}$. The second criterion calculates dissociation energies for all possible fragmentation channels $\text{Si}_{n}\text{O}_{m}\rightarrow\text{Si}_{k}\text{O}_{l}+\text{Si}_{n-k}\text{O}_{m-l}$ with $0\leq k\leq n$ and $0\leq l\leq m$:

[TABLE]

and picks the lowest of them: $E_{\mathrm{diss}}(n,m)=\min_{k,l}\{E_{\mathrm{diss}}(n,m,k,l)\}$.The higher is $E_{\mathrm{diss}}(n,m)$, the more resistant to fragmentation the cluster is. In stable clusters both $\Delta_{\mathrm{min}}(n,m)$ and $E_{\mathrm{diss}}(n,m)$ should be positive, while a negative value of $\Delta_{\mathrm{min}}(n,m)$ or $E_{\mathrm{diss}}(n,m)$ is a sign of instability.

Fig. 3 shows the contour maps of calculated $\Delta_{\mathrm{min}}(n,m)$ and $E_{\mathrm{diss}}(n,m)$ as functions of $n$ and $m$. In both figures, the areas of high stability look like mountain ridges or islands. As expected, silica $(\text{SiO}_{2})_{n}$ clusters are highly stable according to both criteria. Surprisingly, $\text{Si}_{n}\text{O}_{m}$ clusters with $n\sim 2/3m$ exhibit comparable stability. These clusters resemble $(\text{SiO}_{2})_{n}$, but are constructed of Si–O–Si bridges only and have no Si=O double bonds (see Fig. 2a). As seen in Fig. 3a there are also several minor stability islands running along $n\sim 4+2m/3$, $n\sim 6+2m/3$ and $n\sim 10+2m/3$. Such non-stoichiometric compounds are also rather stable according to the second criterion (Fig. 3b). They contain excessive silicon, which segregates as a compact group of Si atoms, attached to the skeleton of Si–O bonds only (see Fig. 2b). Such clusters are of interest due to experiments on growth of long silicon nanowires from gas-phase SiO [24].

Another interesting class of stable compounds is super-oxidized $\text{Si}_{4}\text{O}_{10}$ and $\text{Si}_{n}\text{O}_{2n+1}$ ($n\geq 5$) clusters (see Fig. 2d). The latter have relatively low values of $\Delta_{\mathrm{min}}$, but quite high values (2.5-3 eV) of $E_{\mathrm{diss}}(n,m)$. We also distinguish the \ceSi4O18 cluster which is notably stable by the first criterion. Interestingly, this cluster is similar to recently synthesized \ceP4O18 cluster38, but has free-ending O3 groups instead of closed ones in \ceP4O18 (see Fig. 2d). The important feature of super-oxidized clusters is spin-polarized groups with O-O bonds providing their high reactivity, which may determine the toxicity of silica particles6, 39.

The exploration of $(\text{SiO})_{n}$ clusters is of interest for astrophysics due to the presence of SiO molecules in the circumstellar space and their role in formation of silicates. We note that these clusters do not form a distinct range of $\Delta_{\mathrm{min}}$ or $E_{\mathrm{diss}}$. For this reason they can transform to neighboring, more stable clusters that should be taken into account in interpreting their optical signatures.

Concluding, we have developed a new method for simultaneous prediction of structures of clusters in vast areas of composition. Comparing to currently used methods our approach demonstrated 5-50 times speed-up, allowing for massive ab-initio calculations of nanoclusters at reasonable costs. The availability of such an efficient tool opens the door to wide exploration of trends in chemistry of multi-component nanoclusters, to study cluster features connected with the bulk $x$–$T$ phase diagram, and to the search for new, nonstoichiometric ’islands of stability’, which can be interesting for applications. These prospects for nanomaterials science are supported by our first-principles study of $\text{Si}_{n}\text{O}_{m}$ clusters in a very wide range of compositions. We present the overall picture of stability in these clusters and show numerous ridges and islands of stability, which are very distinct from well-studied silica clusters. We hope, this first attempt gives strong impetus to wide ab initio research in the plethora of important multi-component cluster systems with rich chemistry.

{acknowledgement}

A.R.O. is supported by the Russian Science Foundation (No. 16-13-10459). We also thank Russian Foundation for Basic Research (No 16-02-00612, 18-32-00991 and 17-02-00725). Calculations were performed on the Rurik supercomputer at MIPT, on MVS-10p cluster at the Joint Supercomputer Center (Russian Academy of Sciences, Moscow, Russia) and Lobachevsky cluster at the University of Nizhny Novgorod.