Stability and metastability of clusters in a reactive atmosphere: theoretical evidence for unexpected stoichiometries of MgMOx

By applying a genetic algorithm in a cascade approach of increasing accuracy, we calculate the composition and structure of MgMOx clusters at realistic temperatures and oxygen pressures. The stable and metastable systems are identified by ab initio atomistic thermodynamics. We find that small clusters (M<= 5) are in thermodynamic equilibrium when x>M. The non-stoichiometric clusters exhibit peculiar magnetic behavior, suggesting the possibility of tuning magnetic properties by changing environmental pressure and temperature conditions. Furthermore, we show that density-functional theory (DFT) with a hybrid exchange-correlation (xc) functional is needed for predicting accurate phase diagrams of metal-oxide clusters. Neither a (sophisticated) force field nor DFT with (semi)local xc functionals are sufficient for even a qualitative prediction.

In the search for novel functional materials, atomic (sub)nanometer clusters are widely studied as model systems exhibiting unique size-dependent properties often even qualitatively different from bulk materials. For example, small clusters may exhibit completely new local structures, stoichiometries, electronic and magnetic properties unknown in the bulk materials [1]. Heterogeneous catalysis is just one important example where all mentioned issues are of fundamental importance [2][3][4][5][6][7][8].
The composition and structure of clusters are determined by thermodynamics and kinetics at the relevant temperature (T ) and the nature of the environment. In thermodynamic equilibrium, only structures and compositions that minimize the free energy of the combined gas+cluster system will be stable. Although a system is often not in thermodynamic equilibrium, thermodynamic phase diagrams serve as guidelines and important limits for predicting properties and functions of real materials. In this Letter, we address the issue of stability and metastability using a model system that is relevant for many practical applications: free metal (Mg) clusters in an oxygen atmosphere.
Most of the previous research on clusters focused on properties of stoichiometric (Mg M O x , x = M ) clusters [7,[9][10][11][12][13][14][15][16][17], and only few attempts have been made to study properties of the non-stoichiometric (x = M ) clusters [8,9,15]. However, the decisive issue of stability and metastability of clusters with different compositions at realistic conditions (exchange of atoms with an environment) has not been addressed so far. We consider a wide range of Mg M O x cluster sizes: 1 ≤ M ≤ 15 and x determined by thermal equilibrium with the environment at given temperature T and partial oxygen pressure p O2 . For each stoichiometry, the energy is minimized with respect to both geometry and spin state. Unexpectedly, our results reveal (see Fig. 1) that nonstoichiometric clusters with x > M are more stable at realistic (T, p O2 ) when M < 5, while for bigger clusters there is a competition between stoichiometric and more- conditions. The geometries were optimized with PBE+vdW, the (harmonic) vibrational free energy was evaluated with the same functional, and the electronic energy was calculated using PBE0+vdW. The label of the horizontal axis shows the amount M of Mg atoms, and the amount x of O atoms and spin multiplicity M for the non-stoichiometric MgM Ox cluster thermodynamically most stable at pO 2 = 1 atm and T = 300 K. (the stoichiometric clusters are all singlets, i.e., M = 1). The structure of clusters at selected sizes is also shown (all structures are shown in Suppl. Material) than-stoichiometric composition (x > M ).
For the determination of low-energy structures, we use a genetic algorithm (GA). GA mimics the process of "natural" selection to evolve a pool of atomic structures until the structures that fit best chosen selection criteria are found [18][19][20]. In this work, we search for structures that minimize the DFT total energy within each stoichiometry. Note that GA is not a single method, but a whole class of methods, and must be optimized and validated for each system. Details of our implementation of GA  [39]. For the cases where singlet (M = 1) and triplet (M = 3) spin states are almost degenerate, the two sets of energies are shown (for DFT only, since reaxFF does not describe spin). In these cases, the spin state of MgOx is chosen to be singlet for triplet MgOx+2, and triplet for singlet MgOx+2, to preserve the spin-conservation rule (O2 is always kept in its ground triplet electronic state). The geometry differences of clusters with different spin states are invisible at this scale. are given in the Suppl. Material, and its validation is discussed below.
The free energy as a function of T and p O2 is calculated for the minimum of the potential-energy surface (global minimum, GM) and (energetically) adjacent local energy minima for each stoichiometry using the ab initio atomistic thermodynamics (aiAT) approach [23][24][25][26], recently extended to cluster systems [27,28]. The thermodynamic phase diagram is then constructed by selecting cluster compositions and structures with the lowest free energy as a function of (T , p O2 ).
Reliability of predictions on relative stability of clusters with different structures depends on the accuracy of the underlying potential-energy surface (PES). Here, we find that comparing clusters with different stoichiometry poses an additional challenge, which is not apparent if only one stoichiometry is considered. Fig. 2 shows a comparison of DFT energies of the reaction MgO x + O 2 → MgO x+2 calculated with different xc functionals: generalized gradient approximation PBE [29] and hybrid PBE0 [30], both corrected for the van-der-Waals (vdW) interaction (PBE+vdW and PBE0+vdW), and the highest level currently achievable within the DFT framework, i.e., the renormalized secondorder perturbation theory (rPT2) [31], here applied on PBE orbitals (rPT2@PBE). All DFT calculations are performed with the all-electron FHI-aims package, employing numerically-tabulated atomic orbitals [32]. It can be seen that PBE+vdW strongly overestimates stability of clusters with larger x, resulting in a qualitatively incorrect prediction that O 2 adsorption would be favored over desorption up to a large excess of oxygen. Such behavior is not confirmed by the PBE0+vdW hybrid functional or rPT2@PBE (note that a similar behavior is found for PBE compared to PBE0, i.e., without vdW correction).
Interestingly, for lower O 2 -coverage the difference between PBE+vdW and PBE0+vdW/rPT2@PBE energies of O 2 adsorption on Mg and MgO 2 is small despite the error in the O 2 binding energy (see Fig. 2; the calculated O 2 binding energy is 6.23 eV for PBE+vdW, 5.36 for PBE0+vdW, 4.59 for rPT2@PBE, and the experimental value is 5.21 [33]). This can be explained by error cancellation for the clusters: If an O 2 molecule adsorbs nondissociatively on Mg M O x , the error in the description of Mg M O x+2 will cancel with the error in the description of O 2 when calculating the adsorption energy. Indeed, we find that adsorption of O 2 on Mg and MgO 2 does not lead to breaking of O-O bonds. The difference between PBE+vdW and PBE0+vdW energies of O 2 adsorption on MgO for the triplet case is due to the difference in the description of the singlet state of MgO itself. For clusters with x ≥ 5, correction of the O 2 binding energy error with the experimental value increases the difference between PBE+vdW and PBE0+vdW/rPT2@PBE adsorption energies (see the Suppl. Material). The tendency of PBE (and, as a consequence, of PBE+vdW) to overbind O 2 molecules at high coverage holds also at larger M (see the Suppl. Material).
It has been recently shown [34] that HSE06 functional [35] yields a good description of the level alignment and electron transfer in MgO. Based on comparison for selected clusters, we find that PBE0, which belongs to the same family of functionals, yields formation energies essentially identical to HSE06 (within 0.05 eV for Mg M O x neutral clusters). Taking into account this comparison and the fact that PBE0+vdW results are in general much closer to rPT2@PBE than PBE+vdW, we conclude that reducing the self-interaction error by including Hartree-Fock exchange in DFT is crucial for correct description of the Mg M O x cluster energetics, both within each stoichiometry and in particular when comparing different stoichiometries.
Therefore, the PBE0+vdW functional is used to calculate the energies and evaluate the fitness during GA runs. In order to improve the efficiency of the GA scan, our implementation proceeds in terms of a cascade in which successive steps employ higher levels of theory, with each next level using information obtained at the lower level. Initially, a local optimization of a given structure is performed with lower-level (computationally relatively cheap) numerical settings. At this level, PBE+vdW and "light" numerical settings with basis set "tier 1" were used [32], and forces were converged to better than 10 −3 eV/Å. Next, structures with energies within 2.5 eV from the current GM candidate, are further relaxed using higher-level settings [36]. We use PBE+vdW functional with "tight -tier 2" numerical and basis settings for energy minimization at the higher level, and forces were converged to better than 10 −5 eV/Å. Next, the energy of these further optimized structures are re-evaluated using PBE0+vdW and "tight -tier 2" settings.
A challenging problem of any random-walk-type multidimensional global minimization scheme (including basin hopping and GA) is to guarantee that the lowest-energy structure found by the algorithm is indeed the GM. We address this validation problem by applying replicaexchange molecular dynamics (REMD) [27,37,38], a (computationally very expensive) reference method that performs the canonical sampling of a PES simultaneously at different temperatures, and is exhaustive if the system is ergodic. The validation of GA is performed using a reactive force-field (reaxFF [39,40]) to evaluate energy and forces. While the force-field is by far not accurate enough to predict correct energy differences, as can be seen in Fig. 2 and in Suppl. Material, using it for the validation is a more stringent test for the GA, since the force-field PES is found to have a much more complicated landscape, with many more local minima, than the ab initio PES. With our implementation of GA we find, for the systems here considered, the same GM as found by 1.5-µs-long (cumulative time) REMD [21] runs. An independent evidence of the robustness of our GA scheme is that, for stoichiometric clusters, we could always find the GM reported in literature [3,17].
At given T , p O2 , and M , the stable stoichiometry of a Mg M O x cluster is determined via aiAT, i.e., by minimizing the Gibbs free energy of formation [28]: Here, F Mg M Ox (T ) and F Mg M (T ) are the Helmholtz free energies of the Mg M O x and the pristine Mg M cluster (at their ground state with respect to geometry and spin), respectively, and µ O (T, p O2 ) is the chemical potential of oxygen. F Mg M Ox (T ) and F Mg M (T ) are calculated using DFT information and are expressed as the sum of DFT total energy, DFT vibrational free energy in the harmonic approximation, as well as translational, rotational, symmetry-and spin-degeneracy free-energy contributions. The dependence of µ O on T and p O2 is calculated using the ideal (diatomic) gas approximation with the same DFT functional as for the clusters [28]. Note that µ O (T, p O2 ) is in turn a sum [28] of the same freeenergy contributions as for the free clusters, where the p O2 dependence is captured by the translational free energy term, i.e., k B T ln p O2 +f (T ) (where the second term does not depend on p O2 ). The phase diagram for a particular M is constructed by identifying the lowest-freeenergy structures at each (T , p O2 ). As a representative example, we show in Fig. 3 the phase diagram for M = 4. Phase diagrams based on reaxFF, PBE+vdW, PBE0 (without vdW), and rPT2@PBE can be found in the Suppl. Material. We find that at all DFT levels, the phase diagrams are qualitatively and quantita- tively very similar for T > 200 K and p O2 < 10 −5 atm. For higher pressures and/or lower temperatures, however, PBE+vdW predicts larger x in Mg M O x as thermodynamically more stable, compared to PBE0+vdW and rPT2@PBE. This is consistent with the results shown in Fig. 2, i.e., PBE tends to favor adsorption of a larger number of O 2 molecules. PBE0 and PBE0+vdW yield almost identical phase diagrams. Thus, at the considered cluster sizes the vdW interactions, within the scheme of Ref. [41], do not affect the differences between free energies of formation of competing Mg M O x clusters. ReaxFF-based diagrams were evaluated for the sake of comparison: they are very different from DFT-based ones and no conclusions, even qualitative, can be drawn from the analysis of the reactive force field results.
For creating the phase diagrams, we have approximated the configurational free energy with the harmonic vibrational free energy. If the clusters exhibit fluxional behavior or are melted, as it is the case for some metal clusters, this approximation can become invalid at high or even moderate temperatures [27,42]. For all the Mg M O x clusters described here, we have tested the validity of the harmonic assumption by running ab initio MD simulations for 20 ps at T = 800 K [43]. We found that, while the structures exhibit large vibrations, their topology is never destroyed, and that the initial 0 K structure was reversibly recovered by cooling down the hot sample. At most, some of the O 2 moieties (see below) either freely rotate or roam around the clusters at larger x. In order to quantitatively account for these larger displacements in term of configurational entropy, we have evaluated the excess free energy of selected clusters, by integrating the thermodynamic relation ∂(βF (β))/∂β = U β , where U β is the canonical average of the total energy at temperature T = 1/k B β. The differential equation is numerically integrated by evaluating U β via 4-ps-long (after equilibration) MD trajectories at five increasing temperatures, from T = 10 K (where the reference free energy is safely approximated by the harmonic expression) to T = 800 K. The integration path was checked to be reversible. In this way all anharmonic contributions are included through the canonical sampling of the PES. For the Mg 4 O x , the anharmonic correction enlarges the stability region at higher temperatures of the highest stoichiometry, Mg 4 O 12 , by at most 50 K at p O2 = 1 atm. Note that the (harmonic or beyond) configurational free energy is only one part of the total free energy and that, at higher (T , p O2 ) and stoichiometries x, the free energy is dominated by the term xk B T ln p. For this reason, the plots are extended up to 1 000 K, where weakly bound O 2 moieties may evaporate; nonetheless, at thermodynamic equilibrium the pressure of O 2 molecules ensures that another O 2 molecule will on average replace the evaporated one. Furthermore, we find (see Suppl. Material) that neglecting translational, rotational, (harmonic or beyond) vibrational, and symmetry-and spin degeneracy-related free-energy contributions results only in slight changes in the phase diagrams, comparable to the differences between PBE0+vdW and rPT2@PBE diagrams, for all considered M .
We have constructed phase diagrams for 1 ≤ M ≤ 15. At thermodynamic equilibrium in a wide range around normal conditions, small Mg M O x clusters (M ≤ 5) are found to be preferentially non-stoichiometric (x > M ). For sizes M > 5 we observe a competition between stoichiometric and non-stoichiometric stable structures. In Fig. 1, we also report the relative stability of stoichiometric and non-stoichiometric structures at low (T = 100 K) and high (T = 1 000 K) temperatures. The relative stability of stoichiometric structures increases with temperature.
A further interesting property we find is that Mg M O x non-stoichiometric clusters have in general several near-degenerate (within ∼0.05 eV) electronic states with different spin. The highest multiplicity we predict is M = 7 for the case Mg 4 O 12 , and values of M = 3 and 5 are largely represented. The analysis of the spin density shows that the unpaired electron density is localized on O 2 and O 3 moieties (see examples in the Suppl. Material). The orbitals occupied by the different unpaired electrons in the same cluster have vanishing overlaps due to large separation. The different spin-states for each isomer can be formed and coexist in oxygen atmosphere without violating spin conservation rules, since successive adsorption and desorption of the (triplet) O 2 molecules in the gas phase allows the clusters to reach all the stable spin-multiplicities observed for each cluster. The thermodynamically favored access of oxygen in small clusters is in sharp contrast to bulk pristine MgO, where stoichiometric composition is strongly favored. In all cases where O 2 or O 3 moieties are coordinated to maximum two Mg atoms, these moieties host an unpaired spin. When the number of Mg is increased such that every O atom can be coordinated to three or more Mg, and at the same time the clusters have large HOMO-LUMO gap (i.e., all valence Mg electron are transferred to O x ), then stoichiometric clusters become thermodynamically stabilized.
All the quasi-degenerate spin states are populated at finite temperature, thus all the non-stoichiometric clusters with energetically quasi-degenerate states are paramagnetic. In contrast to non-stoichiometric clusters, we find that stoichiometric structures are always singlet, separated from the higher-multiplicity states by at least 1 eV. Thus, the stoichiometric clusters are diamagnetic. By looking at Fig. 3, we note that in the range of realistically achievable pressures encompassed by the rectangle, non-stoichiometric (paramagnetic) structures are thermodynamically more stable at lower temperatures, while at higher temperatures the stoichiometric (diamagnetic) structures become more stable. The temperature at which this transition occurs is a rapidly varying function of p O2 . We have thus identified a class of systems that undergo an unusual paramagnetic-diamagnetic transition induced by T and p of the reactive atmosphere, where the change in magnetic behavior reflects the change in composition. The transition between paramagnetic and diamagnetic behavior is smooth when environmental conditions are changed smoothly, because different stoichiometries always coexist within few k B T . In fact, the sandcolored areas in Fig. 3 are the regions where the free energies of the competing compositions/structures are within an energy range of 2k B T .
In conclusion, we have presented a theoretical framework for predicting structure and stoichiometry of stable and metastable clusters in thermodynamic equilibrium with a gas atmosphere. An efficient and unbiased scan of the potential energy of the clusters at various compositions is combined with ab initio atomistic thermodynamics, a tool for evaluating the relative free energy of structures by knowing their electronic energy, vibrational frequencies, and structural parameters for the evaluation of translational and rotational entropy. The methodology has been applied to Mg clusters in an oxygen atmosphere. We have shown that small alkaline earth metal clusters form thermodynamically stable non-stoichiometric "nano-oxides", which have been overlooked so far. They are also predicted to have peculiar (T, p)-dependent, magnetic properties.
We acknowledge the cluster of excellence "Unifying Concepts in Catalysis" (UniCat, sponsored by the DFG and administered by the TU Berlin) for financial support.