The stable and metastable, as measured using an all-electron density functional theory approach, stoichiometric clusters of boron, aluminium, gallium, indium and thallium oxide are reported. Initial candidate structures were found using an evolutionary algorithm to search the energy landscape, defined using classical interatomic potentials, for alumina and india followed by data mining or rescaling. Characterization of the refined structures was performed by electronic structure techniques at the hybrid density functional and many-body GW levels of theory. We make accurate predictions of the spectroscopic properties represented by mean ionization potentials of 11.4, 9.9, 9.8, 8.8 and 8.4 eV and electron affinities of 0.05, 1.1, 1.6, 1.9 and 2.5 eV for boria, alumina, gallia, india and thallia, respectively. The changes in the global minima, atomistic and electronic properties with respect to the cluster and cation size are discussed.
Sesquioxides span a wide range of physical and chemical properties and are used as catalysts, lasing and more generally light-emitting materials, elements of complex photoelectric devices, etc. (Catlow et al. 2010). In particular, india is widely used as a transparent conducting oxide in modern optoelectronic devices, ranging from solar cells to flat-panel displays (Edwards et al. 2004), whereas alumina is not only important in many industrial applications (Hart 1990), in particular as a support for catalysts, but also as particles in rocket exhaust where it is possibly responsible for the depletion of stratospheric ozone (Robinson et al. 1997; Fernandez et al. 2003). Questions arise as to the nature of their physico-chemical properties below the bulk regime. Moreover, is it possible to exploit the size–property relationship to further enhance useful properties? The difficulty in structure and property characterization below 5 nm presents an ideal opportunity for theoretical guidance and prediction. In this regime, strong deviations are expected from bulk-like local coordination, and quantum confinement effects will be greatly enhanced. Low-energy atomic structures of (X2O3)n nanoclusters are therefore investigated for n=1–5 and X=B, Al, Ga, In and Tl. There is already a number of published papers on clusters of alumina (Ying et al. 1993; Fernandez et al. 2003; Van Heijnsbergen et al. 2003; Demyk et al. 2004; Karasev et al. 2004; Song et al. 2004; Fernandez et al. 2005; Patzer et al. 2005; Neukermans et al. 2007; Sierka et al. 2007; Feyel et al. 2008; Santambrogio et al. 2008; Sun et al. 2008; Sierka 2010), gallia (Deshpande et al. 2006) and for india (Walsh & Woodley 2010). The atomic structure of small india clusters was predicted by refining plausible low-energy structures that were generated by an evolutionary algorithm (EA) search over the energy landscape. Rather than a direct search of the density functional theory (DFT) energy landscape (Santambrogio et al. 2008; Doll et al. 2010), classical interatomic potentials were used during the initial global searches and DFT during the local refinements. The same efficient approach of switching energy functions is applied here. The initial cluster structures for all five compounds were generated by scaling the low-energy local minimum (LM) structures, for both alumina and india, data mined from our initial global searches and from the literature. Importantly, the relaxation of scaled clusters can lead to new LM configurations not yet found, even when the configurations of the new compounds are expected to be similar to those of the original compound (Sokol et al. 2010; Woodley et al. 2010). The spectra of the electronic levels are then generated for the all-electron DFT refined, lowest energy structures.
(a) Energy of formation
It is assumed that the configuration adopted by the cluster (X2O3)n is that which minimizes its energy of formation. As in previous work (Walsh & Woodley 2010; Woodley et al. 2010), four different energy functions are employed, each assumed to be better than the previous and, importantly, more expensive to compute.
The first energy function is based on the rigid-ion model (RM), where each cation and anion is represented as a point charge of +3|e| and −2|e|, respectively, and the analytical interatomic potential (IP), 2.1is used to model the interaction between ions i and j, which are distance rij apart. The first term is simply the Coulomb contribution (in electronvolts) between point charges qi and qj, and the latter terms are the standard Lennard-Jones and Born–Mayer terms. For the RM, the energy of formation is simply the sum of all two-body interactions defined by equation (2.1). The values for B and C are taken from earlier work where the bulk phases are modelled (Ooms et al. 1988; Walsh et al. 2009). We have included the r−12 term to penalize interaction distances that are less than a typical bond length (A=1.0 eV Å12 for like-charged species and 10.0 eV Å12 otherwise). In order to represent electronic polarization of oxygen, the Shell model (SM; Dick & Overhauser 1958) is employed (our second energy function). In the SM, each point charge from the RM is replaced by two: one represents the nucleus and core electrons, and the other, which is connected via a harmonic spring to the first, the centre of charge for the shell of valence electrons. Typically, the last three terms in equation (2.1) act between the shells, and the first term does not act between two point charges of the same ion. Both models are implemented within the GULP code (Gale & Rohl 2003).
For the third and fourth energy functions, we use an all-electron DFT method with a local numerical orbital basis set, as implemented in the FHI-AIMS code (Blum et al. 2009; Havu et al. 2009). The FHI-AIMS two-tiered basis set is employed for the oxygen and boron atoms, and only the one-tier basis set for the larger cations, as they have a stronger cation nature, which provides an accuracy of −5.39, −2.20, −10.63, −26.19, −24.61 and −15.24 meV per atom for O, B, Al, Ga, In and Tl, respectively. Within the current approach, the eigenstates obtained from the solution of the Kohn–Sham equations are Gaussian smeared with a 0.05 eV width. For our third definition, the PBEsol functional (Perdew et al. 2008) is chosen as it is unbiased and is not too computationally expensive, while scalar-relativistic effects are treated at the scaled zero-order relativistic approximation (ZORA) level (Vanlenthe et al. 1994). In contrast, for our fourth definition, the hybrid PBEsol0 functional and the full ZORA level of theory (Vanlenthe et al. 1994) are employed. The hybrid PBEsol0 functional replaces 25 per cent of the PBEsol electron exchange with exact Hartree–Fock-like exchange (Ernzerhof & Scuseria 1999; Perdew et al. 2008), where the effective exchange–correlation functional becomes 2.2Furthermore, quasi-particle corrections to the PBEsol0 Kohn–Sham eigenvalues () are included through the GW approximation (Hedin 1965; Hybertsen & Louie 1985), which combines the single-particle Green function (G) with the dynamically screened Coulomb interaction (W). The resulting quasi-particle energies () emerge as a first-order perturbation, calculated using Kohn–Sham molecular orbitals (ψn), 2.3which have a direct correspondence with photoemission (N−1 states) and inverse photoemission (N+1 states) or, in the case of the highest occupied and lowest unoccupied molecular orbital (HOMO and LUMO), levels of the nanoclusters of interest, measures of the respective ionization potentials and electron affinities. The GW approximation has had recent successes in describing both the band structure and defect properties of metal-oxide systems (Jiang et al. 2009; Rinke et al. 2009).
Throughout the paper, we will use ‘RM’, ‘SM’, ‘PBEsol’ and ‘PBEsol0’ to indicate which of the four definitions, given above, is employed.
(b) Generating the stable and metastable configurations
The initial aim is to obtain a set of low PBEsol0 energy clusters, including the global minimum (GM), formed from stoichiometric (X2O3)n units.
First, we employ an EA technique (Al-Sunaidi et al. 2008; Woodley et al. 2010), implemented in the GULP package (Woodley et al. 1999), to search for low-energy local minima (LM) on the RM energy landscapes of (Al2O3)n and (In2O3)n. Here, we assume that the set of lowest PBEsol0 energy LM structures can be readily obtained using standard local optimization techniques to refine the low-RM energy geometries obtained from the EA. For each of the two cations and each n, at least 10 runs of the EA is performed, each starting from a different population of 30 cluster structures. The initial cluster structures are composed of 2n cations and 3n oxygen anions placed randomly within a sphere of radius 8.0 Å. Competition between structures within the current population is simulated in order to determine which structures procreate (phenotype crossover and mutation operators used to generate 30 children structures from 30 parents in each EA cycle) and also which survive into the next EA cycle (where the population consists of only the 30 best unique cluster structures). A Lamarckian version of the EA is used, i.e. structures generated in each of the 1000 EA cycles are subject to an immediate local optimization, in which both LM and saddle points on the RM energy surface are generated. The probability of a structure winning (either for procreation or survival) is based purely on a comparison of its relaxed RM energy and that of the other clusters in the current population. More details of our EA implementation to generating stable and metastable low-energy cluster structures can be found elsewhere (Woodley 2004; Woodley et al. 2004, 2010). We also note that the india LM cluster configurations found from searching over the (In2O3)n landscapes have been previously published elsewhere (Walsh & Woodley 2010).
As shown in figure 1, the 10 lowest RM energy candidate structures from each EA run are refined using the SM and standard local optimization techniques, which are implemented within GULP (Gale & Rohl 2003). At the end of this second stage, we have a set of stable and low-energy metastable cluster structures for (Al2O3)n and (In2O3)n, i.e. the approximate configurations to be used as initial geometries in an electronic-structure-based approach. Moreover, as the cluster structures predicted for alumina are found to be similar to those predicted for india, rather than applying the EA to the other compounds of interest (boron, gallium and thallium sesquioxides), we simply rescaled our approximate configurations so that sensible bond lengths (for B–O, Ga–O and Tl–O, respectively) are obtained. Data mining our results (Sokol et al. 2010) for alumina and india is far quicker than employing global optimization techniques for the three other compounds, although we should be aware that important LM may be missed for larger clusters (higher values of n), especially those of the end member (B2O3)n, where the character of the bonding may differ.
The atom positions of the scaled approximate configurations are relaxed so as to minimize the PBEsol energy—again we use standard local optimization techniques, but this time as implemented within the FHI-AIMS code (Blum et al. 2009; Havu et al. 2009). Consider the PBEsol energy landscape. Ideally, to make the comparisons between compounds, for each configuration we would like all five relaxed energy cluster structures. However, if the target configuration corresponds to an LM within a basin, B, that has a small width, then the target LM may be difficult to find, e.g. the PBEsol energy point of the initial approximate configuration may lie outside B and structural relaxations distort the cluster towards a different LM. Moreover, the target configuration may not be stable (or metastable) for all compounds. Where the final structure does not resemble the initial approximate configuration, we: (i) tried rescaling the PBEsol minimized cluster structure of the neighbouring compound (where target LM has been successfully obtained) and (ii) added the new configuration to the set of approximate configurations if not already included (and thus also targeted for the other compounds). Rescaling in (i) also reduces the number of optimization steps.
Finally, to improve the energies and the respective electronic structures, we recalculated these for each PBEsol-relaxed configuration using PBEsol0 (single point) and the GW approximation. From the PBEsol0 energies, the smaller (n=1 or 2) configurations are labelled ‘nR’, where R is the rank, starting from a (the GM) and increasing alphabetically, of the configuration as predicted for the mid-sized cation members, (Ga2O3)n. For larger configurations, n>2, R is composed of either one or four letters that indicate the rank of the cluster for boron oxide or the ranks for the other compounds, i.e. 5abdc indicates that the n=5 cluster is the GM for alumina, ranked second for gallia, ranked fourth for india and ranked third for thallia.
3. Results and discussions
(a) The (X2O3)n clusters for B2O3, Al2O3, Ga2O3, In2O3 and Tl2O3
The three saddle-point configurations shown in figure 2b were generated for (X2O3)n clusters of size n=1. The order of stability (or rank) as measured by PBEsol0 is the same for each cation (figure 2a), i.e. configuration 1a is the GM, 1b the next lowest energy LM and 1c the highest energy configuration. Ignoring alumina, as the (Al2O3)1 configuration 1b is the PBEsol-GM configuration rather than 1a, the energy differences between the three configurations increases when PBEsol0 is employed rather than PBEsol. The energy difference is also greatest between the (B2O3)1 configurations, then (Tl2O3)1, i.e. there is a greater stability of the 1a configuration (compared with 1b and 1c) for the end members. The PBEsol0 energies for these smallest stoichiometric clusters suggest that shorter average bond distances are preferred as opposed to a higher average coordination of the constituent atoms (figure 3).
As reported previously for (In2O3)1 (Walsh & Woodley 2010), cluster 1a relaxes to a linear configuration if an RM is employed, i.e. if polarization effects are ignored. As polarizability of an ion increases with increasing ionic radius, we would expect and do find that the X–O–X and O–X–O angles, α and β, respectively, increase towards 180° as the size of the cation decreases. A linear GM configuration is found for both (B2O3)1 and (Al2O3)1. The latter is in agreement with B3LYP calculations (Sun et al. 2008), which also suggests that larger sized alumina clusters adopt fullerene structures (which is our finding for boria, and not alumina—see below). However, a different DFT study suggests that the alumina cluster should have α≠180° (Song et al. 2004).
The 1b clusters are composed of a tetragonal ring and an additional oxygen atom attached to one of the cations (forming an additional stick). The tetragonal ring is the (XY)2 GM structure for many XY compounds, for example, (ZnO)2 (Al-Sunaidi et al. 2008), (ZnS)2 (Hamad et al. 2009) and (MgO)2 (Roberts & Johnston 2001), whereas the additional stick is the trivial GM for (XY)1. The observation that a (X2O3)1 cluster can be constructed from a stable (XO)2 cluster provides a possible method for generating (X2O3)n clusters from (XO)2n clusters, i.e. one can construct (X2O3)n clusters from (ZnO)2n global minima configurations found from a previous study (Al-Sunaidi et al. 2008) and decorate the outer cations with the remaining n oxygen anions (with a maximum of one additional oxygen anion for each cation) to create potentially low-energy (X2O3)n configurations. Below, such clusters will be noted.
Of the 1b configurations, (In2O3)1 has a tetragonal ring that is most like a square: O–In–O and In–O–In angles are α=92.2° and β=89.7°. Crudely, these angles gradually change with the size of the cation: for the smallest cation size of boron, α=134.0° and β=68.9°, and for the largest cation of thallium, α=82.0° and β=96.0°. A similar trend is found for configuration 1c: as the rate of change of interatomic distance X–X is more rapid than X–O (as clearly seen in figure 3a), α decreases and β increases with cation size. In contrast, the O–X–O angles, γ, formed by the additional oxygen anion, are approximately 136° for all compounds (which would be 120° if there was no second cation). The sum of α and β is not expected to equal 180° due to the presence of this additional oxygen anion. There are some discrepancies from these trends if we look more closely at the values of α or β of 1b for alumina and gallia. Although the respective angles were similar, which relates to the similar sizes of aluminium and gallium cations, α(Al) is in fact less than α(Ga), and β(Al) is predicted to be greater than β(Ga).
The HOMO and LUMO energies are shown for the (X2O3)1 configuration in figure 2. For each compound, the between these energies, i.e. EGAP=ELUMO−EHOMO, are greatest for the GM configuration 1a. With the exception of B2O3 configuration 1b, EGAP as measured using PBEsol is quite small for the metastable configurations and so smearing occupation of eigenlevels was problematic (ideally, smearing should be switched off). However, such problems do not occur when the GW approach is employed as EGAP increases for all configurations, with HOMO becoming more stable and ELUMO increasing. Note that HOMO is degenerate for the linear arrangement of configuration 1a.
(b) The (X2O3)2 clusters for B2O3, Al2O3, Ga2O3, In2O3 and Tl2O3
The lowest energy structures (one per unique configuration) for the n=2 (X2O3)n clusters are shown in figure 4. We have labelled them according to their rank as found using PBEsol0 for gallia. In figure 5a, 3.1the PBEsol0 binding energies relative to the energy of n non-interacting 1a GM configurations, are shown for each configuration i. Missing data points correspond to configurations for which we were unable to find a corresponding local stationary point on the PBEsol energy hyper-surface, i.e. on local optimization, the configuration (for example, scaled from that found for a neighbouring compound) no longer resembled the initial structure. As already discussed in §2b, the energy basin containing the LM may either have a very small or zero width, thus making it hard or impossible to find. For some initial configurations, where structural changes are probably owing to the cation size, we show more than one relaxed cluster structure with the same label: the O–X–O–X chain in configuration 2f is bent if composed of cations larger than Al3+, whereas there is less bonding (fewer sticks shown in figure 4) in configurations 2g, 2h and 2l for (B2O3)2 and likewise for 2m composed of cations smaller than In3+. Although configurations 2d and 2k have the same connectivity, we labelled them separately as both configurations were easily found for the larger cations—(Tl2O3)2 prefers the planar arrangement of 2k, unlike the other compounds. In the (B2O3)2 example of configuration 2l, the cluster has actually fragmented into two (B2O3)1 clusters.
We now return to our earlier proposal, in §3a, of constructing LM cluster structures by adding n oxygen atoms to LM (ZnO)2n clusters. Two very different low-energy configurations for 1–1 compounds are found for compounds like ZnO: an octagonal ring and a cuboid (Al-Sunaidi et al. 2008). Adding two oxygen atoms to these, we form configurations 2m and 2n, respectively. Although LM, if we ignore , both were found to be quite high in energy.
For each compound, the energies in figure 5a are arranged using the rank found for the corresponding configuration for (Ga2O3)2, thus, the energy always increases for the gallium oxide clusters. The indium oxide clusters have similar relative energies, between −6 and −2 eV, to those of gallium oxide, likewise, but between −2 and 0 eV, for boron oxide and thallium oxide. The relative energies for aluminium oxide are significantly lower, which follows the trend found for the linear 1a configuration with respect to configuration 1b (see figure 2, where the energy difference is smallest for Al2O3, similar for Ga2O3 and In2O3, and much larger for the end members B2O3 and Tl2O3). Clearly, the linear configuration is less favourable for aluminium oxide, and as a consequence, there is a greater driving force for aluminium oxide to form larger clusters (moreover, the planar configurations 2k and 2m were either not found or much higher in energy). Although we have not investigated the minimum energy path to forming an n=2 configuration from two isolated 1a configurations (and thus predicting energy barriers, if they exist), we can imagine two linear clusters coalescing to form configuration 2l, which, for alumina, is approximately 2 eV lower in energy than the planar configurations and still approximately 2 eV higher than most other n=2 configurations and more than 4 eV higher than the GM configuration.
The energy differences, , of the n=2 configurations, where is the average of for the set of (X2O3)2 configurations, are plotted as a function of cation (from the smallest, boron, to the largest, thallium) in figure 5b. Examining this figure, significant changes in the rank of configurations can be clearly seen (typically associated with a steep gradient). For example, configurations 2a, 2b and 2i have many crossovers from boron to aluminium. More importantly, we can see that 2c and 2b are the GM configurations for (B2O3)2 and (Al2O3)2, respectively, and 2a otherwise. Although much higher in energy for most compounds, the planar configuration 2m is the second lowest energy configuration for (Tl2O3)2. A previous DFT study for alumina clusters (Song et al. 2004) considered a number of configurations including configurations 2a (their GM) and 2c, we assume that they missed configuration 2b.
The ionic radii of the (six-fold coordinated) cations are 0.41, 0.68, 0.76, 0.94 and 1.03 Å. As Al3+ and Ga3+ have similar sizes (difference of less than 0.1 Å), we might expect both the configurations and PBEsol0 ranking to be more similar than what our results predict. For example, alumina configurations 2i and 2j have a much better ranking than that predicted for gallia as they are lower in energy than configurations 2d–h. The GM configuration for alumina, 2b, is also different—2a is the GM configuration for X=Ga, In and Tl. Thus, it again appears that the more compact clusters (with higher coordination numbers) are more stable for alumina. In a similar way, with a difference of ionic radii of less than 0.1 Å, we can also find some similarities of configurations and ranking for (In2O3)2 and (Tl2O3)2; for example, a different 2m configuration is found for other compounds, and E2n is significantly higher than the energy of other configurations reported here (whereas is lower than ).
As already mentioned above, the n=2 GM configuration changes with cation size. Bond lengths and angles of the three configurations found are shown in figure 6. As found for the n=1 clusters, the bond lengths increase with cation size, while the tetragonal rings become more square (angles nearer 90°) and angle α in configuration 2a widens (although always remains below 180°). The bond angles within the 2b, however, remain constant at approximately 114° (120° would be the perfect angle for a three-coordinated site). Interestingly, the search over the RM energy landscape for alumina generated an additional LM cluster with the same connectivity as 2b, but flattened with a bond angle, β, of 107.5° and 142.7° (93.9° and 134.5° in the SM) that had D2d point symmetry rather than Td. This second lowest RM energy configuration relaxes to the GM configuration 2b when the PBEsol energy is minimized.
(c) Atomistic structure and interatomic potential energetics for (Al2O3)n and (In2O3)n (n=3−5)
The search over the rigid ion (Al2O3)3 energy landscape yields LM structures that are essentially the same as those previously reported (Walsh & Woodley 2010) for (In2O3)3. Using the labels in figure 7 to distinguish the different configurations, 3aaaa was found to be the GM for (Al2O3)3, then, 3bbbb, 3cddh, 3aaaa, 3eccc, 3dejk, 3ij-i, 3fihe, 3–fd, 3*ggif, 3klll, 3-feg, 3hhg- and then 3jkkj. The last three structures were data mined from the configurations generated for india (so there is a greater chance of finding low-energy LM configurations not shown in figure 8), and an asterisk is used to indicate that the structure changes symmetry upon relaxation: C2v to Cs, with atoms in the upper half (as seen in figure 8) moving to the left so that the initially four-coordinated central oxygen atom is only three coordinated. The rank, as measured using the RM, for india are similar: 3aaaa, 3cddh, 3bbbb, 3eccc, 3dejk, 3ij-i, 3fihe, 3*ggif, 3–fd, 3-feg, 3hhg-, 3klll and then 3jkkj. Again, the symmetry of 3ggif is lowered when relaxed, and the last two configurations were data mined from the LM of alumina. The three lowest energy LM structures for both compounds resemble a tea cosy with symmetry point groups C1, C2v and Cs. The C1 tea-cosy structure is the rigid ion GM for both compounds, whereas the most symmetric, open tea-cosy structure is ranked second for india and third for alumina. For alumina, we found more than one C1 tea-cosy LM: ranked first and fourth. Similarly, as reported in §3d, when relaxing these configurations on the PBEsol energy landscape for all five compounds, more than one C1 tea-cosy configuration is found. As the bonding (indicated in figure 8 using sticks) and symmetry are the same, below we only report the energy of the most stable configuration for the C1 tea cosy. Further similarities are found between the LM ranking of the two compounds. Considering the order of the india configurations ranked 4–8, a D3h bubble and four Cs structures (Walsh & Woodley 2010), we only need to reverse the order of the last two to obtain the ranking for alumina. This change in rank suggests that tetragonal faces are not penalized as much for alumina. In fact, a rhombic dodecahedron, with Oh symmetry, an oxygen atom at its centre and not previously reported for india, is the next lowest energy LM configuration (3klll) found on the (Al2O3)3 rigid-ion energy landscape. In this configuration, the six aluminium atoms form an octahedron about the inner oxygen atom, and each aluminium atom is also coordinated to four of the eight outer oxygen atoms, forming 12 edge-sharing Al2O2 tetragonal faces. For india, there are at least two more configurations with a better rank than 3klll, including the more planar, less compact structures, 3-fee and 3hhg-.
The n=4 and 5 lowest energy stationary point structures found by searching the RM landscapes for (Al2O3)n and (In2O3)n after refinement on the PBEsol landscape are also shown in figure 7 together with a few high-energy structures already reported in the literature (Sun et al. 2008; Walsh & Woodley 2010). The initial order, or rank, of these clusters for (Al2O3)4, starting from the GM, is 4caaa, 4egeg, 4gikl, 4ddbb, 4jhde, 4kjcd, 4bchi, 4fegf, 4lkih, 4i-jf, 4abfc and 4hf-k, where the last five configurations were initially generated using data mining, whereas on the (In2O3)4 IP landscape, it is 4caaa, 4egeg, 4gikl, 4ddbb, 4bchi, 4jhde, 4kjcd, 4i-jf, 4fegf, 4hf-k, 4abfc and 4lkih, with two of the last three configurations obtained from data mining our relaxed PBEsol (B2O3)4 configurations. As already noted, the configurations in figure 7 were obtained after minimizing the PBEsol energy for one of the five compounds: alumina (cations coloured purple in the online version), boron oxide (pink) and indium oxide (grey). Note that we have chosen to show one particular configuration twice: configuration 4a has the same connectivity (bonds) as 4hf-k, but is shown rotated by approximately 90° about the vertical axis. As well as a rescaling of bond lengths and bond angles, which is more significant if the cations are changed, there are a few more noticeable differences between the IP and PBEsol-minimized configurations. Ignoring structural changes that involve bond breaking, the largest difference is found for configurations 4i-jf and 4hf-k: a cation–oxygen–cation bond angle points inward (towards the centre of the cluster) for PBEsol, but points outwards for IPs (configuration 4hf-k becomes 4a).
For the largest clusters shown in figure 7, the order, or rank, again starting from the GM, for (Al2O3)5 on the IP-landscape is 5ddca, 5vvpr, 5opm-, 5lqoo, 5ieed, 5cjhn(C2), 5hlnk, 5egd-, 5mhki, 5wwrp, 5pkje, 5gnij, 5cjhn(C2v), 5jffc, 5fcaa, 5fcab, 5kigf, 5qmlg, 5rtvs, 5bbs-, 5tstq, 5aabh, 5suw-, 5c, 5uru-, 5a and finally 5- - -m. The first 10 configurations were generated from a direct search over the IP-landscape for (Al2O3)5. And for (In2O3)5, it is 5ddca, 5opm-, 5ieed, 5cjhn(C2), 5egd-, 5lqoo, 5mhki, 5hlnk, 5vvpr, 5cjhn(C2v), 5gnij, 5jffc, 5qmlg, 5wwrp, 5fcaa, 5rtvs, 5kigf, 5aabh, 5bbs-, 5pkje, 5suw-, 5tstq, 5noql, 5uru- and then 5- - -m, with the first eight configurations generated by a direct search over the IP-landscape for (In2O3)5. Of these configurations, large structural distortions occur during the relaxations of 5gnij, 5qmlg and 5- - -m (after switching from PBEsol back to IP). Configuration 5cjhn represents two distinct stationary points on the IP landscapes: the higher symmetry configuration shown in figure 7 has two central four-coordinated oxygen atoms that are bonded to the lowest eight cations that are part of the four outer hexagonal faces. The cluster with C2 symmetry can be obtained from the C2v cluster by displacing atoms vertically so that the coordination of the internal oxygen atom is lowered to two, and two of the four edge-sharing hexagonal faces become pairs of edge-sharing tetragonal faces. On the alumina and india IP-landscapes, the higher symmetry configuration is actually a saddle point, and is higher in energy than the C2 configuration. Configurations 5a and 5- - -m are lower symmetry versions of configurations 5uru- and 5bbs-. The latter, 5bbs-, which appears (as viewed in figure 7) to be similar to 5suw-, is in fact an open cone pointing into the page, whereas configuration 5suw- is a bubble.
Although the IP GM configuration is the same for both alumina and india, i.e. 5ddaa, there are some noticeable differences in the ranking of the other LM, especially for 5vvpr and 5wwrp. Both configurations have one oxygen atom within a cigar-shaped bubble composed of 14 tetragonal and four hexagonal faces. For 5wwrp, the internal oxygen atom bridges two cations, whereas in 5vvrp, five more tetragons are formed (by the four new X–O bonds). Again, we suggest that tetragons are penalized more for clusters containing the larger indium atoms, as 5vvrp is ranked second for alumina and ninth for india, and 5wwrp also falls four places from 10 to 14.
Upon using the SM or an electronic-based method, i.e. including polarizability of the oxygen atoms, as reported previously for ZnO clusters (Al-Sunaidi et al. 2008), two-coordinated anions at the surface move outwards (sharper Al–O–Al bond angles). Also, comparing O–Al–O bond angles for an aluminium atom at the surface that is coordinated to three oxygen atoms that are also at the surface (see, for example, 5fcaa, 5noql and 5tstq), a more planar arrangement (sum of bond angles nearer to 360°) is found in the PBEsol LM configurations than in the IP LM; typical examples include 354.7° and 350.8° (PBEsol0), which are significantly higher than 324.3° and 315.8° (IPs). This trend is of course only valid if no bonds are broken or created.
(d) Density functional theory energetics for (X2O3)n (n=1–5 and X=B, Al, Ga, In and Tl)
As already discussed above, LM configurations on the PBEsol energy landscape do not always resemble the initial structure (before relaxation) that are data mined from the LM configurations on the IP energy landscapes. This discrepancy occurred most frequently for (B2O3)n, where the larger clusters relaxed typically into very different configurations, which are more open as pairs of edge-sharing faces become a single face (bond, or shared edge, is broken). As such, the labels used in figure 7 to distinguish the configurations refer to either the rank found only for (B2O3)n or the ranks found for (Al2O3)n, (Ga2O3)n, (In2O3)n and (Tl2O3)n. Of the remaining four, it is expected that the thallia clusters would produce the greatest structural distortions of the original dataset of structures. Depending on how configuration 3ggif is initially rescaled, the cluster composed of the larger cations will relax to 3ggif or 3–fd. For boria, the uppermost tetragon of 3ggif rotates by 90° about the axis that goes through the two uppermost boron atoms (thus breaking two bonds, while keeping C2v symmetry).
Configuration 3dejk has the point symmetry of D3h and the lowest PBEsol0 energy for (B2O3)3. As the size of the cation is increased, the rank of this configuration falls: fourth for (Al2O3)3, fifth for (Ga2O3)3, tenth for (In2O3)3 and eleventh for (Tl2O3)3. Moreover, the similar ranks match the similarity of the cation ionic radii of Al3+ and Ga3+, and In3+ and Tl3+. A falling (or improving) rank with cation size is more the exception than the rule: 4bchj and 4gikl (4ddbb, 4lkih, 5ddca, 5kigf, 5pkje, 5qmlj and 5wwrp) are the other obvious examples in figure 7. Configurations 3dejk, 4bchj and 4gikl contain large (greater than hexagonal) faces, however, 4hf-k does not have a fall in rank between (Al2O3)4 and (Ga2O3)4, as might be expected, as it is very similar to 3dejk. One interesting (X2O3)3 cluster is that of the dodecahedron that, compared with the other clusters in figure 7, is ranked last on the PBEsol energy landscapes (compared with fourth and second from last on the alumina and india IP landscape); with the shortest metal-inner oxygen bond length, the (B2O3)3 dodecahedron has more of a cubic appearance (outer oxygen atoms are at the corners of a cube).
Configuration 4a has the point symmetry of C2v and the lowest PBEsol0 energy for (B2O3)4. The effect of increasing the size of the cation was investigated; upon relaxation (PBEsol minimized after rescaling), configuration 4hf-k is obtained (an oxygen atom moves within the cluster—see earlier discussion), with the X–O–X angle increasing from Al to Tl, although always remaining smaller than that for B. A similar change, movement of an oxygen atom from outside inwards, is found for 4i-jf, in which the bond angle Al–O–Al (pointing outwards) is 158°, whereas for In–O–In and Tl–O–Tl (pointing inwards), it is 150° and 152°—the direction of this angle was not found for (Ga2O3)4 and is, perhaps, the cause of the difficulty of finding this configuration for gallia.
Configuration 5a with the lowest PBEsol0 energy for boria can be obtained from configuration 5uru-, whereas configuration 5b is essentially 5bbs-. Apart from three tetrahedral boron sites in the latter, the boron atoms in both structures are coordinated with three oxygen atoms. The only other configurations shown in figure 7 are 5c and 5d, which were readily found for boria. For the largest cation, data mining produced the lowest energy configuration twice. Namely, configuration 5ddca has one bond less than 5fcaa for alumina (Al–O distance of 3.61 Åcompared with 2.00 Å, where a typical bond is 1.77 Å), whereas for thallia, the two configurations are identical (Tl–O distance of 2.68 Å compared with a typical bond of 2.03 Å).
The PBEsol0 energies of the configurations shown in figure 7 are plotted twice in figure 8. Typically, the energy of formation is reported relative to the n=1 GM configuration, as in figure 8a. From the predicted structures, we can conclude that the boria clusters are generally unlike the configurations adopted by the other compounds, particularly as n increases. In contrast, figure 8a indicates a similarity between gallia and india, and boria and thallia, with alumina clusters being much more stable (the thick line has the same profile for each n—a minimum for Al). As the energies of clusters with similar-sized cations are likely to be similar, this observation does not meet expectations—in fact, the trend highlighted by the thick line in figure 8a is merely an artefact of the relative stability of the n=1 cluster with respect to typical larger clusters. Therefore, we plotted the relative PBEsol0 energies again, but this time using the symmetrical 2b configuration (point symmetry Td), where the cations and anions have a coordination of three and two, respectively. The mirror image of the earlier profile is now found for n=1. For n>1, the change in the relative energy is significantly different for the boria GM configurations, as expected. The interesting result is what happens to the other compounds: initially (n=1–3), stability of the GM increases with cation size, but by a decreasing magnitude with increasing n. For n=4, similar GM stability (with respect to its smaller sized cluster) is predicted for alumina and gallia, and similarly for india and thallia, whereas for n=5, a greater stability is found for alumina (and by a smaller magnitude india).
The high end of the highest density of the light plus symbols in each column in figure 8 does not always provide an estimate of the transition between the highest energy configuration produced using global optimization techniques and the lowest energy configuration that were obtained by data mining. Although the worst ranked n=5 configuration for larger cations, 5rtvs, was generated by data mining the literature, for smaller cations, the worst ranked configuration reported here, 5wwrp, was generated by the global search over the alumina IP landscape (although ranked tenth).
Ignoring clusters of boria, the same GM PBEsol0 configuration is predicted for n=3, namely the C1 tea cosy 3aaaa. This structure has also been predicted for alumina and india on the respective IP landscapes. The tea-cosy structure is also reported for (Al2O3)3 in a study using ab initio total energy linear combination of atomic orbitals (LCAO) calculations (Fernandez et al. 2003). For n=4, configuration 4caaa, with point symmetry D3d, is the PBEsol0 GM for gallia, india and thallia, and can be considered as a fragment cut from corundum (bulk phase). For alumina, this configuration is the third lowest LM, after 4abfc and 4bchi—both of which rank well for the clusters containing similar-sized cations, viz. Ga. The corundum fragment was previously reported as a possible GM for alumina (Song et al. 2004) and has been used when modelling water on (Al2O3)n using ab initio total energy LCAO calculations (Fernandez et al. 2003); however, configuration 4abfc is also reported as the alumina GM structure by other groups using EAs to search the DFT (B3LYP) energy landscape (Sierka et al. 2007; Sierka 2010). For n=5, cations of a similar size have the same GM PBEsol0 configuration: alumina and gallia adopt 5aabh, whereas for india and thallia, it is 5fcaa—a maximum change in rank that is caused by a cation size of seven (a–h). In figure 8, more than one black cross is shown in a particular column if the same GM (X2O3)n configuration is not predicted for Al, Ga, In and Tl. The distance, or energy difference, between the dark plus symbols, therefore, indicates when a change in configuration can help reduce the strain caused by a change in cation size. The largest change is found for (In2O3)4 between configurations 4caaa and 4abfc, where the change of rank for india is five.
Energy differences between the two lowest energy configurations for each compound and cluster size are shown in figure 9. The tea-cosy structures of n=3 (first four configurations in figure 8) have similar energies and result in the minimum for the clusters containing the larger four cations; at room temperature these tea-cosy structures are likely to be part of one ergodic region on the free-energy landscape, thus we have also included the energy difference between the GM and the lowest energy non-tea-cosy n=3 structure. Now the energy differences in figure 9 show that the relative stability of each GM configuration with respect to clusters of the same size generally decreases with n. Again, it is interesting that there is not more similarities between the results shown for alumina and gallia; ignoring the boria clusters, the energy difference for (Ga2O3)n is similar to (In2O3)1, (Tl2O3)2, (In2O3)3 and (Tl2O3)5.
The complexity of the energy landscapes will typically increase with the number of variables or atoms. Stability is, of course, not just dependent upon the energetic barriers (Schön et al. 2003). Kinetic, or entropic barriers will become increasingly important as the number of atoms, or size of the clusters, increases.
(e) The electronic structure for (X2O3)n (n=1–5 and X=B, Al, Ga, In and Tl)
Analysis of their electronic states (single-particle PBEsol0) reveals that each of the identified stable and metastable (X2O3)n structures has finite HOMO–LUMO separations, with the highest energy occupied states being dominated by O 2p and the lowest lying empty states consisting largely of X Ys and Yp orbitals for (X,Y)=(Al,2), (Ga,3), (In,4), (Tl,5). The density functional calculations are therefore consistent with the formal charge states of X3+ (Ys0Yp0) and O2− (2p6), which explain the excellent correspondence between the classical and quantum energy landscapes for these larger cations. As expected, the eigenvalues resulting from the semi-local PBEsol functional typically result in a small HOMO–LUMO gap, whereas the non-local PBEsol0 functional, which counters the self-interaction error and the energy derivative discontinuity present in the pure DFT functional (Sham & Schlüter 1983), opens the gap through both a significant lowering of the HOMO level and a raising of the LUMO level. Previously, we reported that the application of the GW quasi-particle correction for india GM clusters maintains the HOMO position of the PBE0 functional (to within 0.3 eV) (Walsh & Woodley 2010). Here, for the PBEsol0 functional, we typically obtained a similar initial HOMO position, but now the application of the GW quasi-particle correction lowers the HOMO position; for n=3–5 by approximately 1.2 eV. This movement is dependent upon the cation; it is larger for alumina, approximately 2 eV, and smaller for thallia, approximately 1 eV. For all compounds, the LUMO further raises towards the vacuum level by approximately 1–2 eV.
The GW spectra near the HOMO (first peak below −5 eV) and LUMO (first peak above −5 eV) levels of the ground state (GM) clusters for all compounds and sizes investigated here are plotted in figure 10. The eigenvalue for HOMO, or quasi-particle ionization energy, for boria is −12.5 eV for n=1, −10.0 eV for n=2, −11.5 eV for n=3 and 4 and −11.3 eV for n=5. More interestingly, three of these clusters have a positive LUMO eigenvalue suggesting that they would not accept an additional electron from the vacuum, and the electron affinities for the other two are small: 0.2 eV for n=3 and 0.04 eV for n=5. The quasi-particle ionization energies (electron affinities) in electronvolts for alumina GM clusters n=1–5 are −9.9 (1.6), −8.9 (0.8), −10.3 (1.0), −10.3 (1.1) and −10.2 (1.1). Similar values of the ionization energies are obtained for the first four gallia GM clusters: −10.0 (1.4), −8.9 (1.4), −10.0 (1.7), −9.9 (1.8) and −10.0 (1.9). The quasi-particle ionization energies (electron affinities) in electronvolts for india GM clusters n=1–5 are: −9.1 (1.7), −7.8 (1.5), −8.9 (1.2), −8.9 (2.3) and −9.1 (2.6). These are more similar to that obtained for thallia: −9.1 (1.5), −7.1 (1.9), −8.5 (2.9), −8.6 (2.9) and −8.6 (3.1). The mean values of the ionization potentials therefore decrease with cation size (11.4, 9.9, 9.8, 8.8, then 8.4 eV), whereas the electron affinities increase (0.05, 1.1, 1.6, 1.9 and 2.5 eV). The energy difference between HOMO and LUMO levels is, for a particular size n, typically largest for boria and smallest for india and thallia. Moreover, for each n, the spectra for india and thallia are similar, which is not surprising as they adopt the same GM configurations, whereas boria adopts very different GM clusters for n>1, and the spectra for boria clusters are less like those of the other compounds.
Since the same GM configuration is not found for all compounds, in figure 11 the spectra for two low-energy clusters, configurations 2b and 4caaa, are plotted for each compound. As expected, the spectra for the same configuration are similar. The position of the LUMO level is typically constant, whereas the position of HOMO and the density of states for the clusters containing larger sized cations typically increase. Now consider the non-GM cluster configurations that are generated from the GM configurations after changing cations, e.g. 4caaa for (X2O3)4 is one such LM if X=Al. The spectra of these non-GM clusters can be added to those of the GM clusters. If the spectra are also weighted with the Boltzmann factor, , where Na is the number of atoms in the GM or LM cluster and kBT=0.025 eV, then we obtain the room temperature spectra in columns two and four in figure 11 for n=2 and 4, respectively. Ideally, all LM clusters that would add a significant contribution to the room temperature spectra of the electronic energy levels should be included. Aligned in figure 12 are the spectra for 10 LM (Al2O3)4 clusters, including the GM, the resulting room temperature spectra for (Al2O3)4, and the observed spectra for bulk alumina. The room temperature spectra for (Al2O3)4 are remarkably similar to those observed for bulk alumina.
The stable and metastable, as measured using an all-electron DFT approach, stoichiometric (X2O3)n clusters of boron, aluminium, gallium, indium and thallium oxide have been reported. Initial candidate structures were found using an EA to search the energy landscape, defined using classical interatomic potentials, for alumina and india followed by data mining, or rescaling (Sokol et al. 2010). Electronic characterizations of the refined structures were performed at the hybrid density functional and many-body GW levels to obtain accurate predictions of the spectroscopic properties. For n>1, (B2O3)n adopts different GM configurations than that predicted for the other compounds. With similar ionic radii, alumina and gallia were expected to adopt similar GM configurations. However, for n=2 and 4, the alumina GM configurations were different from that of gallia, which adopted the same GM configuration predicted for india and thallia. Comparisons of cluster structures indicate that X2O2 tetragons are not as heavily penalized for the smaller Al3+ cation. Also, for the n=1 clusters, boria and alumina adopt a collinear structure, whereas the greater polarizability of the larger cations results in the buckling of this structure. The relative energy of a cluster is typically either reported with respect to the smallest sized (n=1) GM cluster or to the bulk. Here, we show that the use of the former is misleading when making comparisons of relative energies of different compounds. Finally, the electronic spectra near the HOMO and LUMO levels were reported for two example configurations and the GM configurations for all five compounds. The mean values of the ionization potentials for the GM clusters decreased with cation size, from boria to thallia, whereas it increased for the average electron affinities. Furthermore, a Boltzmann (room temperature)-weighted superposition of spectra for (Al2O3)4 LM was found to be similar to the experimental data taken from the X-ray photoelectron spectrum and the ultraviolet photoelectron spectrum for bulk alumina.
The author is grateful for valuable input from Alexey Sokol, Aron Walsh, Stephen Shevlin and Richard Catlow, and for funding from UK’s HPC Materials Chemistry Consortium, which is itself funded by EPSRC (EP/F067496). In the completion of this work, the author made use of molecular visualizing software from Accelrys and the UCL Legion High Performance Computing Facility, and associated support services; and the facilities of HECToR (approx. 1 million AUs), the UK’s national high-performance computing service. The latter is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc and NAG Ltd, and funded by the Office of Science and Technology through EPSRC’s High-End Computing Programme.
One contribution of 16 to a Special feature ‘High-performance computing in the chemistry and physics of materials’.
- Received January 6, 2011.
- Accepted February 9, 2011.
- This journal is © 2011 The Royal Society