We review recent developments and applications of computational modelling techniques in the field of materials for energy technologies including hydrogen production and storage, energy storage and conversion, and light absorption and emission. In addition, we present new work on an Sn2TiO4 photocatalyst containing an Sn(II) lone pair, new interatomic potential models for SrTiO3 and GaN, an exploration of defects in the kesterite/stannite-structured solar cell absorber Cu2ZnSnS4, and report details of the incorporation of hydrogen into Ag2O and Cu2O. Special attention is paid to the modelling of nanostructured systems, including ceria (CeO2, mixed CexOy and Ce2O3) and group 13 sesquioxides. We consider applications based on both interatomic potential and electronic structure methodologies; and we illustrate the increasingly quantitative and predictive nature of modelling in this field.
Computational methods have a long history of applications in the study of materials used in energy technologies; and as shown by several articles in this issue, they continue to play a substantial role and one, moreover, which is likely to grow in the future. They are applied to the study of materials used in both energy generation—including nuclear, battery, fuel-cell and solar energy systems—and storage—including gas storage, most notably of hydrogen. This article will concentrate on modelling materials at the atomic level, where their role is to provide information on the key atomic level structures, processes and parameters that control the behaviour of the material in its applications in either generation or storage. In particular, substantial contributions have been made by modelling methods to our knowledge of the following.
— Defect structures and properties including defect formation and migration, which control transport properties and may determine both spectroscopic and thermodynamic processes. Early examples include studies of the fission reactor fuel, uranium dioxide (Catlow 1977, 1978, 1987; Jackson et al. 1986; see also other articles in the special issue of Journal of the Chemical Society, Faraday Transactions 2, vol. 83, issue 7, 1987), which established the basic defect structure of the material including the energies for defect generation and for both cation and anion migration. Simple calculations also showed that there was a relatively low-energy inter-cation charge transfer process (Macinnes & Catlow 1980), which was proposed to be responsible for the high temperature specific heat excess in the material—an important phenomenon in possible accident scenarios as discussed in greater detail by Stoneham (2010) in this issue. Defect processes also play a crucial role in battery and fuel cell materials; closely related is the modelling of intercalation which is of major importance in battery technologies. And defects processes are of central importance in radiation damage as discussed in the article of Duffy (2010).
— Surface properties including modelling of surface structures, defects and impurities and molecular adsorption and reaction, which are again of major importance in the performance of materials in energy generating or storage devices. In this context, it should be noted that catalytic processes play a central role in several energy technologies.
— Spectroscopic properties which are clearly vital in solar energy materials and where defects may often provide the light-absorbing centres.
— Molecular sorption and diffusion, relating to materials used in gas storage, of both hydrogen and CO2.
Modelling methods may also be used to calculate a wide range of relevant properties, including elastic and dielectric constants. The methods may be used to predict new materials as well as to model properties of existing systems; and indeed, the predictive capacity of atomistic modelling continues to grow.
Modelling of energy materials makes complete use of the extensive range of techniques available in contemporary condensed matter computational science. Methods based on interatomic potentials—minimization, Monte Carlo and molecular dynamics—all have a substantial and continuing role in the field; and there have been very extensive efforts over the last 30 years in developing and refining potential models particularly for oxide materials used in energy technologies. The last 15 years has seen a very rapid expansion in the use of electronic structure methods, based on both density functional theory (DFT) and Hartree Fock (HF) techniques, with the former being particularly prominent. Detailed discussions of these methods are available elsewhere (Martin 2004; Dovesi et al. 2005; Yip 2005; Jensen 2006; Catlow et al. 2010). All will continue to be needed in investigating the complex and varied problems posed by energy materials.
In the sections which follow, we will present illustrative applications of both achievements of and problems posed by modelling of energy materials. We will mainly highlight the work of the present authors and their collaborators and we note that other examples of successful applications can be found elsewhere in this issue, especially in the articles of Islam (2010), Duffy (2010) and Stoneham (2010).
Our review will focus first on the materials used for clean energy generation and storage, specifically, for hydrogen production by water-splitting. The challenges involved in hydrogen storage necessitate research into new classes of materials reviewed next. The latter naturally leads to the key issues of alternative energy storage and utilization in fuel cells and batteries. Materials used in direct conversion between solar power or light and electricity will be explored thereafter. We conclude with a general account of the search for stable nanostructures, as nanostructured materials are of key importance in contemporary energy generation and storage technologies.
(a) Materials for hydrogen production by water-splitting
In this section, we first discuss recent developments in a range of ternary materials with potential applications in water-splitting. We then focus on one particular anode system, SrTiO3, which illustrates some of the challenges in undertaking accurate modelling studies of perovskite-structured oxides. This section concludes with an account of the interaction of hydrogen with cuprite-structured Cu2O and Ag2O, which are prospective cathode materials.
(i) Photocatalysts for water-splitting
Our current utilization of solar energy is limited, even though it is the most abundant and clean energy source available. While direct conversion of visible light to electricity in photovoltaic devices is a desirable process, the generation of a renewable chemical feedstock, such as hydrogen from water powered only by sunlight, would also be highly beneficial. The origins of photoelectrochemistry can be traced back to Becquerel’s discovery of the photoelectric effect in 1839, where illumination of metal halide solutions between two platinum electrodes produced an electrical current (Becquerel 1839). The feasibility of using a photo-electrochemical process on a large scale requires the discovery of efficient and long-lifetime catalysts, with sustainable components, for which metal oxides would be the ideal candidates. Unfortunately, most metal oxides have optical band gaps that lie outside the visible range (Eg>3 eV) and thus, even if they are catalytically active for a given reaction, they cannot make efficient use of the solar spectrum. Based on an average terrestrial solar spectral irradiance, figure 1, a material with a wide band gap such as TiO2 (Eg=3.2 eV) is limited to 1 per cent solar to hydrogen conversion efficiency for the water-splitting process, while for a smaller band gap material such as Fe2O3 (Eg=2.2 eV), the theoretical efficiency is increased to 15 per cent (Rocheleau & Miller 1997). To remain competitive with solar cell powered electrolysis of water, conversion efficiencies of at least 10 per cent will be required.
Considering the water-splitting process in more detail, an n-type photocatalyst would be the active anode, directly facilitating the oxidation of H2O to O2, as illustrated in figure 2. Absorption of light greater than the band gap of the material will generate an electron–hole pair; the electric field present at the semiconductor–electrolyte interface will facilitate charge carrier separation, with the hole (h+) being driven towards the interface:2.1while a p-type photocatalyst (or a working electrode such as Pt or WC) would act as the cathode, forwarding the reduction of H+ to H2 using the photogenerated electron (e−):2.2completing the overall water-splitting reaction:2.3The thermodynamics of these reactions will depend on the relative position of the semiconductor band edges (hole and electron energies) with respect to the (H+/H2) and (H2O/O2) redox levels. The valence band of the semiconductor must be below the (H2O/O2) oxidation level (5.7 V below the vacuum level) to make reaction (2.1) thermodynamically favourable, while the conduction band position must be above the (H+/H2) level (4.5 V below the vacuum level) for reaction (2.2) to occur spontaneously. If both conditions are not fulfilled simultaneously in a single material, an external bias voltage may be applied, or a tandem device could be constructed using two or more active materials to generate charge carriers with sufficient energy to complete the overall reaction.
In some cases, undesirable competitive processes may also occur, e.g. in Cu2O, owing to its high valence band position, which is close to the (H2O/O2) oxidation level (Xu & Schoonen 2000), Cu oxidation to produce CuO is known to take place,2.4as well as the dissolution of Cu ions:2.5
Furthermore, even if the correct band alignment exists for the redox processes to take place, the kinetics of electron transfer on the semiconductor surface may require the presence of a co-catalyst such as Pt, NiO or RuO2 to accelerate the reaction rates. Excellent reviews on these topics are available (Bak et al. 2002; Kudo & Miseki 2009).
Here, we will discuss the progress in the understanding of photoactive metal oxides, taking our simulation work on the newly synthesized, low band gap metal oxide Sn2TiO4 and the prototypical photocatalyst SrTiO3 as examples, as well as discussing future pathways for computational material research in this area.
Band gap engineering: the case of Sn2TiO4
Titanium dioxide is one of the most widely studied metal oxide systems owing to its extensive applications, ranging from heterogeneous catalysis to dye-sensitized solar cells (Linsebigler et al. 1995). TiO2 is a wide band gap n-type semiconductor and is generally, under ambient conditions, in one of three polymorphs: rutile, anatase and brookite. Tetravalent Ti possesses a formal 3d0 configuration, with the unoccupied 3d states forming the conduction band. For visible light photocatalysis, one wants to sensitize TiO2 to visible light, but for the water-splitting reaction, this must be through raising the valence band, as opposed to lowering the conduction band, because the conduction band already lies at the limit of the (H2/H+) reduction potential, while the valence band is situated well below the (H2O/O2) oxidation potential (Xu & Schoonen 2000).
Contraction of the electronic band gap by raising the valence band could be achieved through modification of the anion site, by introducing ions less electronegative than O with lower binding energy p states, such as S or N, thereby forming the corresponding stoichiometric oxysulphide or oxynitride systems. Modification of the cation sites to incorporate dn transition metals can result in increased visible light absorption, but may also suffer from self-trapped polaronic charge carriers (Walsh et al. 2008, 2009a; Kleiman-Shwarsctein et al. 2009). Alternatively, one could consider cation engineering through the incorporation of species with occupied low binding energy d10 states, e.g. Cu(I) and Ag(I). A more unconventional approach is to introduce species with a cationic lone pair, i.e. group 13–15 cations with filled s bands. Forming multi-component oxide systems with both a d0 cation (e.g. Ti(IV), V(V), W(VI)) and an s2 cation (e.g. Sn(II), Bi(III), Sb(III)) has resulted in a family of stable, visible light absorbing photocatalysts including BiVO4 (Walsh et al. 2009c), Bi2WO6 (Fu et al. 2005), SnWO4 (Cho et al. 2009) and SnNb2O6 (Hosogi et al. 2008). For these systems, the top of the valence band is influenced by the s2 cation, while the bottom of the conduction band is determined by the d0 cation.
Engineering the band gap will also have direct consequences on the defect properties of a material, and hence its conductivity. In general, for an oxide material to support excess electrons, there should be a reducible cation present to provide charge compensation, e.g. Ti(IV) can readily be reduced to Ti(III), making TiO2 n-type, while to support electron deficiency, an oxidizable cation should be present, e.g. Cu(I) can be easily oxidized to Cu(II), contributing to the p-type behaviour of Cu2O (Scanlon et al. 2009). Combining oxidizable and reducible cations in the same system should enable a material to support both electron and hole carriers, and thus facilitate separation of the excitons created on photon absorption.
Following these arguments, the combination of Sn(II) 5s2 and Ti(IV) 3d0 cations in a single material would be desirable. Tin monoxide (SnO) features a small band gap (Eg=2 eV) and has been gathering recent attention owing to its good p-type conduction properties (Ogo et al. 2008). Divalent Sn has an occupied 5s2 state around 7 eV below the top of the valence band, which hybridizes with O 2p and nominally empty Sn 5p states to produce a local distortion in the Sn electron distribution and coordination environment, resulting in the layered litharge structure (Walsh & Watson 2004). The oxide formed from tetravalent tin (SnO2) adopts the rutile structure, similar to TiO2, and exhibits a wider band gap (Eg=3.6 eV) and intrinsic n-type conductivity (Freeman & Catlow 1990; Godinho et al. 2009).
Early attempts to combine Sn and Ti in a ternary oxide system resulted in anatase or rutile solid solutions containing the Sn(IV) species, i.e. without the presence of the occupied Sn 5s2 orbitals (Yamane et al. 1992; Hirata 2000). However, Kumada et al. recently reported the successful synthesis of Sn2TiO4 single crystals from TiO2 and SnCl2 precursors (Kumada et al. 2009). The resulting crystal structure is isostructural to Pb3O4 (red lead) and is shown in figure 3. Pb3O4 is a mixed oxidation state lead oxide, with a crystal structure (P42/mbc) containing two cation sites characteristic of the Pb atoms found in both litharge Pb(II)O (the 8h Wyckoff position with distorted fourfold coordination to oxygen) and rutile Pb(IV)O2 (the 4d position within distorted oxygen octahedra; Wells 1984). The structure is generally considered as a member of the spinel class of oxide minerals, although its crystal structure is unusual and is only adopted by a handful of other oxides including Sb2CoO4, Pb2SnO4 and Sb2ZnO4. There are four formula units per unit cell (28 atoms) and a notable hollow in the centre of the crystal that has been considered to be stabilized by van der Waals interactions between the Pb lone pair electrons (Dickens 1965). The synthesis of Sn2TiO4 in the Pb3O4 structure is therefore not unexpected, given that it simultaneously satisfies the coordination preferences of both the Ti(IV) and Sn(II) species.
The equilibrium structural parameters of Sn2TiO4, calculated using DFT with the PBE functional and a plane wave basis set, employing the VASP code (Kresse & Furthmüller 1996a,b), are within 1.5 per cent of experimental reports. At this level of theory, the band gap is underestimated at 1.16 eV, but using the HSE06 functional, which includes 25 per cent of non-local electron exchange, a separation of 2.09 eV is predicted between the valence and conduction bands. For both Sn and Ti, the calculated Bader partial charges of +1.36 and +2.28 are lower than the formal oxidation states, but are of similar magnitude to those calculated for the separated binary oxides litharge SnO (+1.25) and rutile TiO2 (+2.29); it should be remembered that absolute values of partial charges derived from delocalized charge densities are far from conclusive (Catlow & Stoneham 1983; Jansen & Wedig 2008).
While the valence band is dominated by O 2p states, the presence of Sn 5s is noted at −7 eV, and to a lesser extent at the top of the valence band (figure 4). The mixing of Sn 5p is also significant, consistent with previous analysis for SnO (Walsh & Watson 2004) and PbO (Walsh & Watson 2005). Contributions from Ti dominate the conduction band due to the low energy of the Ti 3d states, similar to TiO2. The band structure exhibits some degree of dispersion at the top of the valence band; the highest states at the Γ and X points are doubly degenerate and composed of Sn 5s and O 2p character. The electron density node between Sn and O, as shown in figure 5, is indicative of the antibonding nature of the interaction at the top of the valence band, which gives rise to the asymmetric electron density (lone pairs) on Sn, which is directed into the hollow of the crystal lattice. Long-range dispersive interactions involving the lone pair electrons are not treated fully at the PBE–DFT level, and while these may further lower the cohesive energy, they will not affect the pertinent structural and electronic characteristics. The conduction band minimum is centred at the Γ point, and composed of predominately Ti d character with some hybridization with O p. The hole and electron effective masses along the Γ to Z line are found to be 1.5me and 0.9me, respectively, which are quite reasonable for polar materials and indicate a better potential for n-type electrical conductivity.
Sn2TiO4 demonstrates that band gap reductions can be achieved through combining s2 and d0 cations, which will be active in facilitating hole and electron transport, respectively. While such low band gaps can be found for transition metal oxides containing dn cations such as Co, the small polaron-dominated conductivity severely limits applications that rely on the generation of a robust photocurrent (Walsh et al. 2009a). Dispersion in both the valence and conduction band edges is found to be moderate for Sn2TiO4, indicating the potential for bipolar conductivity; however, it is unclear whether this material will be intrinsically n-type (as TiO2) or p-type (as SnO), or whether the low band gap will result in good visible light photocatalytic behaviour. These are issues that should be addressed in future experimental and theoretical studies. The present work, however, shows the value of such calculations in guiding the choice and optimization of materials.
(ii) Polymorphism and redox properties of anode material SrTiO3
Now we turn to a more extensively studied multifunctional metal oxide system, SrTiO3. Applications of SrTiO3 as a photocatalyst are matched by interest in its ferroelectric and semiconducting properties. The conduction band level of SrTiO3 is situated above the (H+/H2) redox level (Xu & Schoonen 2000), making it a suitable material for water splitting; however, similar to TiO2, at 3.2 eV its band gap is too large to absorb a significant fraction of visible light, and much effort has gone into engineering the optical absorption properties through chemical doping. For example, the modification of oxygen vacancy concentrations through aliovalent cation substitutions (Takata & Domen 2009) and band gap modification through the incorporation of noble metals (Konta et al. 2004) have both shown promise for increasing the photocatalytic activity under visible light conditions.
SrTiO3 has the perovskite crystal structure. Rotation and tilting of the perovskite building blocks gives rise to a number of polymorphs, which are accessible under specific thermodynamic conditions: the cubic perovskite structure (space group 221, ) as shown in figure 6, a tetragonal expansion containing an octahedral rotation about the c-axis (space group 140, I4/mcm), and a recently proposed orthorhombic modification involving octahedral rotations about other axes (tilts), within space group 74, Pbnm (Smith in preparation). The high-symmetry cubic phase is stable at room temperature and has a transition point at approximately 110 K to the tetragonal phase. To explain the experimentally observed dielectric behaviour of SrTiO3 in the low-temperature regime, a transition to an orthorhombic phase has been proposed to take place at 65 K (Lytle 1964).
Although there are limited experimental data on the orthorhombic structure (Lytle 1964), the tetragonal phase at low temperatures is believed to display rotational distortions of the octahedra of approximately 2° (Tsuda & Tanaka 1995). It was assumed that the rotational and tilting distortions in the orthorhombic structure would be similar (much lower, for example, than in barium titanate). Existing shell model pair potentials predict significantly higher rotation angles of the order of 8–10°. To study the low-temperature orthorhombic strontium titanate we have developed three new sets of pair potentials using a previous set of potentials as a starting point (Akhtar et al. 1995); results of parametrization are summarized in tables 1 and 2.
Potentials a and b were developed to produce two different sets of rotation angles but are otherwise very similar. Potential a produces an optimized orthorhombic structure with octahedral rotation angles of 5.29° and 2.07° and for potential set b 4.88° and 1.59°.
Potential set c was developed to reproduce the structure of potential set a, but with much smaller opposing charges on the titanium core and shell. The large charges and strong spring constant were believed to be responsible for convergence problems in a series of more sophisticated hybrid QM/MM calculations.
Here, we compare properties of SrTiO3 calculated using new potentials with those available in the literature (Akhtar et al. 1995; Crawford & Jacobs 1999). Much of the experimental data and modelling work has been on the room temperature cubic phase, and so for the best comparison of results, the cubic phase has also been used in the present study. All calculations were performed using the experimental lattice constant of 3.9051 Å, with symmetry constrained relaxation of the internal positions.
As shown in table 3, the three new potentials reproduced the experimental elastic and dielectric properties reasonably well, with the exception of the static dielectric constant for potential set c, which indicates an instability arising from a doubly degenerate imaginary frequency of 49i cm−1. One should consider that these potentials were developed to model the orthorhombic unit cell and can be expected to predict instability of the cubic case on the potential energy hypersurface. In this comparison, we follow the previous work (Akhtar et al. 1995; Crawford & Jacobs 1999) where imposing the experimentally observed lattice parameters allows one to model high-temperature phases using potential energy-minimization techniques. By relaxing the lattice parameters, we could study the symmetry breaking which would result in the stabilization of the alternative low-temperature phases. In this context both potentials a and b prove to be sufficiently robust, while potential c fails and should be retained for simulations including explicit account of the thermal motion of the lattice (e.g. using quasi-harmonic approximation or in molecular dynamics simulations), at least in the room-temperature regime; it is also, as noted, more robust when used in QM/MM calculations. From the three models, potential a gives the best overall agreement with the experiment. The Crawford potentials produce similar results, with less accurate elastic constants, a good static dielectric constant and slightly underestimated high-frequency dielectric constant. Akhtar potentials produce more accurate elastic properties and high-frequency dielectric constant than the current work, but underestimate the static dielectric constant by approximately 25 per cent. Finally, of the five sets considered potential set c best describes the orthorhombic phase, as could be expected. Simulations of the low-temperature phase will be discussed in detail elsewhere.
Next, we focus on defect properties of strontium titanate, and in particular, oxygen vacancies, which proved to be the active sites involved in a wide range of physical and chemical processes that control the utilization of this material, including its potential application in photocatalysis.
Formation energies of intrinsic defects
It has been demonstrated that Schottky defect formation is the dominant mechanism of intrinsic ionic disorder in perovskites (Lewis & Catlow 1986; Akhtar et al. 1995) owing to the close-packed nature of the structure. The full Schottky quintet can be written as2.6with the associated defect reaction energy defined with respect to isolated ion removal and the addition of one formula unit to the bulk lattice (E[SrTiO3]):2.7Owing to the high energetic cost of creating a highly charged Ti vacancy, partial Schottky formation in the SrO sub-lattice can become a competitive process:2.8where the defect reaction is defined with respect to the standard state of SrO (the relevant energy is calculated using experimental room-temperature lattice parameters of its rocksalt phase):2.9The association energy (for strontium and oxygen) is the energy of interaction between a strontium and oxygen vacancies being taken from infinite separation:2.10
All defect energies have been calculated here using the Mott–Littleton approach as implemented in the GULP package (Gale & Rohl 2003). In this approach, a site of interest within the crystal, which usually contains a point defect and includes all atoms up to a certain cut-off radius (making region 1) is treated most accurately with explicit relaxation while the effect of the defect on the remainder of the crystal is considered approximately using linear response techniques (up to a further cut-off radius that defines region 2a with an atomic resolution, and further—up to infinity—region 2b which is considered effectively as a continuous medium). The defect calculations used region sizes of 10 Å and 25 Å for regions 1 and 2a, respectively.
The energies of formation for full and partial Schottky defects, as well as the component elementary defects, are shown in table 4 for the three potential models. The partial Schottky defects have been calculated using lattice energies of rocksalt SrO and rutile TiO2 with lattice constants set at their experimental values, as the potentials we discuss are still too simplistic to be transferable between structures with differently coordinated ions.
Our calculations show a good quantitative match to the previously calculated lattice energies and single defect energies. In the case of the full and partial Schottky defect formation energies, our values are in between the previously reported values, though the partial Schottky defect energies are about 0.7 eV short of the experimentally reported values. The association energy for the partial Schottky defect formation is lower than expected from experimental reports (close to zero), but is consistent with the high dielectric constant that dampens the electrostatic attraction between charged defects.
Oxygen and strontium ion migration
Oxygen vacancies are undoubtedly present in SrTiO3 and their diffusion is a matter of significant interest. The diffusion barrier is determined by the lowest energy pathway between two equivalent oxygen positions. To a first approximation, the saddle point lies in the (110) direction, in a linear path equidistant between two oxygen sites (Akhtar et al. 1995; Crawford & Jacobs 1999); however, some curvature in the trajectory is expected. To model this process, we follow the ideal linear path and perform a transition state search close to the centre to find the true saddle point. Strontium diffusion is also of interest whereas, due to the high energy cost associated with its high-charge state, titanium diffusion does not take place to any significant level.
Experimental values for the oxygen migration barrier range from 0.65 to 1.27 eV (Sakaguchi et al. 2000). However, most modern techniques suggest that an O activation barrier of approximately 1 eV is correct. As can be seen from table 5, a range of values for the activation barrier can be obtained using the different potential parametrizations, but potential set c most closely matches the activation barrier energy. The values presented by Akhtar and Crawford reproduce better the lower experimental values. The calculations presented here showed that the oxygen transition point is slightly off the centre.
All of the strontium migration barriers presented here fall directly in the range set in previous simulations and offer a marked improvement in relation to the experiment. The corresponding transition points using all the potentials presented in this work occur directly in the middle between the two vacant strontium sites as illustrated in figure 7.
To summarize, we have developed three new sets of pair potentials and shown that they characterize the cubic system very well in comparison with the experiment and improve upon the previously available pair potentials. Preliminary work on defect formation and transport using these models produces consistent results, and acts as a starting point for performing higher level QM/MM calculations, which are currently in progress. Overall, the work reported here emphasizes the care that is needed in developing potential models for perovskite-structured oxides, and the valuable role that such models can play in characterizing key defect parameters.
(iii) Hydrogen in cathode material Cu2O
As outlined above, hydrogen evolution in a photoelectrochemical cell occurs on the cathode, for which Cu2O is considered to be one of the best candidates due its band gap of 2.2 eV spanning the optical range and its intrinsic p-type conductivity. Cu2O has also gathered interest as a photovoltaic absorber in low cost solar cells (Sears & Fortin 1984). In fact, cuprous oxide was the first material recognized as a semiconductor and also the first material where a Wannier–Mott exciton with a fully resolved optical spectrum has been observed. It is also of interest to compare its properties with the isostructural silver(I) oxide Ag2O, which has a lower band gap of 1.3 eV (Tjeng et al. 1990). Both the band structure and defects in these materials generated much interest in the last 30 years (Nolan 2008) and the most recent defect investigations have been performed at a hybrid exchange and correlation level (Scanlon et al. 2009). As hydrogen evolves by trapping of photoexcited electrons on the surface proton sites, atomic hydrogen can either associate to form gaseous H2—a desired product—or it can diffuse into the cathode material, possibly altering its properties. In our studies we aimed to establish the fate of the atomic hydrogen inside the bulk cuprite lattice.
To find the most stable configurations adopted by hydrogen in these materials, we have carried out a series of density functional Γ-point calculations (using a PBE exchange and correlation density functional along with a double numerical with polarization basis set as implemented in the code DMol3; Delley 1990, 2000) on a 2×2×2 supercell (including Z=16 formula units) of the cuprite lattice. In these simulations, hydrogen was initially placed at all interstitial sites with distinct coordination within the lattice, and then all internal coordinates were allowed to relax. Analysis of the resultant configurations helped us to identify three low-energy configurations of particular interest: (i) a tetrahedral bridging quasi-atomic interstitial; (ii) a bond-centre hydroxyl interstitial; and (iii) an inverted tripodal hydroxyl. These are shown in figure 8.
The charge state of hydrogen—neutral quasi-atom in configuration (i) and protons in configurations (ii) and (iii)—was determined by the analysis of the spin density, which shows for the latter two states a complete transfer of the electron from the hydrogen to the first and second neighbouring copper sites, where it is smeared over cation valence p (for Cu) and dz2 (for Ag) like states. In fact, the spin density from the hydrogen atom in configuration (i) also reveals a strong hybridization between H 1s and nearest neighbour Cu p or Ag dz2 states and to a smaller extent with O sp3 states. Moreover, in Ag2O the electron transfer to the nearest neighbours is nearly complete. Relative energies of the three stable configurations for the two oxides are given in table 6. We used atomic hydrogen as a reference point (rather than molecular), as it is the form in which it appears in the material both in photoelectrochemical cell and muon experiments discussed below.
Using the raw calculated energies, one could conclude that both materials have a strong preference for the inverted tripodal configuration, in which the proton breaks one of the metal oxide bonds by pulling an oxygen ion from its tetrahedral lattice site and leaving the metal ion near its original position with only one bond remaining to the lattice oxygen. At the chosen level of theory, however, we do not expect to obtain very accurate energetics but rather a semi-quantitative estimate of the energy ordering and, perhaps more reliably, structural information. Nevertheless, we can use our knowledge of the inherent problems of the approximations used to produce better estimates of the hydrogen incorporation energies, also given in table 6. Indeed, as the band gap is strongly underestimated by our model (at 0.7 eV for Cu2O and only 0.1 eV for Ag2O), the electron lost by hydrogen in this configuration to the conduction band is overstabilized by a similar amount of ca 1.5 eV for Cu2O and 1.2 eV for Ag2O, which we use to correct the raw calculated data. Thus, we ascertain that the bond-centre configuration in both materials is least likely to exist. We also observe a similar effect for the most stable (as calculated), inverted tripodal hydroxyl configuration and to a lesser extent the quasi-atomic configuration, where the overstabilization is roughly 50 per cent of that of the proton states. From our resulting estimates of the energy of hydrogen incorporation given in table 6, we can conclude that the atomic tetrahedral bridging interstitial sites provide the most stable sites for hydrogen in Cu2O, whereas in Ag2O hydrogen is significantly more stable as the inverted tripodal hydroxyl. The energetic order of the metastable states is also inverted between the two materials. The much narrower gap in Ag2O contributes to the observed ordering of the configurations and is responsible for the different behaviour of hydrogen in this material compared to Cu2O: it has been reported that it is a deep acceptor in Cu2O and a shallow donor in Ag2O (Cox et al. 2004, 2006). Further work is clearly needed using more sophisticated and accurate methods, as in addition to the band gap underestimation, electron and hole localization are severely penalized by GGA-type exchange-correlation density functionals, which can be traced back to the problem of electron self-interaction. This problem is alleviated by hybrid functionals, which have, however, become available in plane-wave codes only recently and still remain prohibitively computationally expensive in most cases. Moreover, to obtain quantitatively accurate description of charged defects or even charge-neutral species, but with a strong dipole, we should employ large supercells thus aggravating the task even further. Therefore, appropriate embedded cluster hybrid QM/MM approaches may be necessary, as for example in our studies of ZnO described below. Nevertheless, we consider that the predicted stable configurations are valid.
Thus, our calculations, in agreement with experiment, suggest that atomic hydrogen remains a very stable species when trapped at tetrahedral bridging interstitial sites, which explains why this material can efficiently convert protons (trapped as surface hydroxyls) to gaseous hydrogen without competing back reactions.
Once the problem of hydrogen production is overcome, we need to address the issue of storing hydrogen in an efficient and sustainable manner as discussed below.
(b) Hydrogen storage materials
In order to reduce anthropogenic CO2 emissions, it is necessary to reduce fossil fuel consumption. One way in which this may be done is to make a switch to cleaner fuels, of which hydrogen is the favoured choice. However, as discussed further by Züttel et al. (2010), hydrogen is difficult to store safely at low cost, light weight and low volume. Thus storage in solid materials that have fast H2 release and adsorption kinetics at 85°C (a typical operating temperature) is desired (Schlapbach & Züttel 2001). No single material meets all of these requirements. For example, magnesium hydride (MgH2) has long been considered as a promising material for hydrogen storage owing to its substantial thermal stability, low cost, the reversibility of the reaction and significant gravimetric hydrogen storage capacity (7.6 wt%), but suffers from a number of drawbacks, including high desorption temperatures (between 300 and 400°C), hydrogen outgassing at a low pressure of 10 Pa (or 10−4 bar) at ambient, and relatively slow absorption/desorption kinetics. Lithium borohydride (LiBH4) reversibly desorbs hydrogen in two steps, following melting at 560 K,releasing 13.1 wt% hydrogen, but as this is a multi-step process with the associated slow kinetics, and as this occurs at heightened temperatures, unmodified LiBH4 is not ideal for applications. Therefore, research in this field is active and ongoing. Simulation can help one to understand and modify the properties of hydrogen storage materials. In particular, the strength of H2 binding is of importance for hydrogen storage applications. This is defined aswhere ETot(XH2) is the total energy of the original system, ETot(X−H2) is the total energy of the original system after two hydrogen atoms are removed and ETot(H2) is the total energy of a H2 molecule in vacuum, chosen as the reference state. In our calculations all geometries are relaxed. The more negative EBind(H2), the stronger hydrogen is bound. In order to meet the US Department of Energy targets for hydrogen storage applications, the binding energy (sometimes referred to as the removal energy or equivalently the absorption energy) should be −0.7 eV/H2≤EBind(H2)≤−0.2 eV/H2 (Shevlin & Guo 2009). The binding energy is a typical property of interest for the hydrogen storage community as it is relatively easy to calculate. Free energies of H2 release are a more accurate measurement of hydrogen affinity, but for large nanostructures it is prohibitively expensive to calculate vibration modes, and the electronic binding energy is a suitable measure for hydrogen affinity.
In this section, we shall illustrate the role of simulations in aiding hydrogen storage research by considering two systems: (i) ideal and defective boron nitride and carbon and (ii) transition metal-doped boron nitride and carbon. All calculations presented are DFT calculations as implemented in the VASP code. Unless otherwise stated, the GGA exchange correlation functional was used.
(i) Carbon and boron nitride systems
Many hydrogen storage materials suffer from slow adsorption and desorption kinetics relating to the need to bring two hydrogen atoms together to form the neutral H2 molecule. Therefore, systems that bind H2 via a physisorption mechanism, i.e. layered compounds such as carbon or boron nitride, are from this point of view preferred. Both boron nitride and carbon share similar structural motifs. Because of the similarity of the B−N and the C−C bonds, boron nitride forms several phases that are isomorphous to carbon systems: i.e. layered hexagonal boron nitride (h-BN) which is similar to graphite and the cubic zinc-blende phase which is similar to diamond. The hexagonal structure can also be rolled up into nanotubes that are structurally similar to the carbon nanotubes (CNTs), but are electrically and energetically different. The partially ionic nature of the B−N bond opens a band gap in the hexagonal boron nitride sheet, whereas the purely covalent nature of the C−C bond allows the hexagonal graphite sheet to be a semimetal. The nature of the B−N bond also affects the structure of closed-shell BN clusters; they do not form fullerene cages as is the case for the carbon clusters as it is not energetically favourable for B−B and N−N bonds to form, as will be discussed in greater detail later.
Initial interest in these materials for hydrogen storage applications was spurred by the early results of Dillon et al. who observed a hydrogen uptake of between 5 and 10 wt% for single-walled CNTs at low temperature (133 K) and very low pressures (0.04 MPa; Dillon et al. 1997). Later experiments however have not reproduced these results. More recent results indicate that CNTs store (at 300 K) approximately 0.1 wt% at atmospheric pressure (Lueking & Yang 2004) rising to 1 wt% at gas overpressure of 1 MPa (Shen et al. 2004).
Boron nitride nanotubes possess an experimental hydrogen storage capacity of 1.8–2.6 wt% at a pressure of 10 MPa at 300 K, but also have an enhanced gravimetric storage capacity at higher temperatures when compared with carbon systems (Ma et al. 2002). The addition of metal nanoparticles such as platinum causes the creation of voids and other point defects in the tubes that boost the hydrogen storage capacity to approximately 4.2 wt% (Tang et al. 2002). Clearly, for both types of material it is important to understand the interaction of hydrogen with ideal sheets and native defects in order to improve their performance. In particular, the bond between B and N atoms is heteropolar in nature; this variation in atomic charge through the material may increase the binding of H2 to the substrate, as the hydrogen molecule is weakly polarizable and possesses a small quadrupole moment. Additionally, the presence of two elements with differing electronegativities may have implications with regard to the splitting of H2 to form atomic hydrogen. We therefore performed calculations to determine the binding of H2 to planar hexagonal boron nitride (h-BN) systems. We used a single h-BN sheet with a 32-atom orthorhombic unit cell of dimensions 8.7×10×18 Å (with a defect density of 3.125%). The optimized in-plane parameters of the primitive surface cell are very close to the experimental value of 2.505 Å (Paszkowicz et al. 2002). As boron nitride is a wide band gap semiconductor, we sample only the Γ-point.
In the lowest energy site for H2 physisorption, the molecule is situated 3.35 Å above the middle of the hexagon, with a binding energy of −0.04 eV/H2 (a similar value is obtained for graphene). In addition, we found that the lowest energy site for dissociative chemisorption is when the H atoms reside on neighbouring B and N atoms, causing them both to adopt partial sp3 hybridization (elongating the B−N bonds and causing the B and N atoms to pop out of the surface by 0.756 Å and 0.579 Å, respectively) with a B−H bond of length 1.234 Å, the N−H bond of length 1.034 Å and the H−H distance of 2.034 Å. This structure is 1.92 eV higher in energy than the contrasting H2 physisorption structure, and thus chemisorption is strongly endothermic by 1.88 eV/H2.
The minimum energy pathway for H2 dissociation on a perfect h-BN sheet is shown in figure 9a. There is a significant barrier of height approximately 0.5 eV preventing chemisorbed hydrogen association to H2; thermal fluctuations will therefore not cause spontaneous dissociation. However, the barrier for H2 to 2H dissociation is much larger (ca 2.4 eV, the same as for graphene). Therefore, it is unlikely that at room temperature H2 will chemisorb onto undefective boron nitride sheets. The transition state has a H−H bond length of 1.166 Å above the B atom, with B−H distances of 1.372 Å and 1.462 Å. One H atom is also bound to the N atom with an N−H distance of 1.197 Å.
By comparison with, for example, graphene, there are two possible single vacancy structures for BN: (where a single B atom is removed) and (where a single N atom is removed). Furthermore, we can consider a double-atom vacancy or bound Schottky pair which we denote , whereby single B and N atoms are removed, which has the benefit of not strongly affecting the charge carrier concentrations and thus preserving the p- or n-type nature of the material. The energies of defect formation have been calculated with respect to standard states (solid boron for , N2 for , and h-BN for ) using standard thermodynamical data (Shevlin et al. 2008). We find that the vacancy formation energies are 10.29 eV , 4.87 eV and 7.30 eV . The interaction of these defects with both molecular and atomic hydrogen species was calculated, and the results are presented in table 7. The reaction energy (Ereac) is defined as the difference in total energy between atomic adsorption and molecular adsorption, the more negative the more energetically favourable is the latter. The single- and double-atom vacancies strongly favour atomic adsorption, but can also form metastable structures, whereby the H2 is physisorbed above the vacancy.
The minimum energy path for H2 dissociation on the three vacancies was determined (see figure 9b–d). The vacancy with the lowest barrier towards molecular dissociation (0.47 eV) is the defect, while the defect with the largest barrier is the defect. This correlates with the length of the H2 molecule at the transition state. For the defect the H2 is stretched by 8.5 per cent (0.813 Å), for the defect the H2 is stretched by 20.8 per cent (0.906 Å) and for the defect the H2 is stretched by 50.9 per cent (1.132 Å). It is clear that exposed nitrogen atoms (created via a boron vacancy) are the most reactive with H2.
(ii) Transition metal doped boron nitride and carbon systems
It is clear from our calculations that unmodified boron nitride and carbon systems do not bind hydrogen strongly enough to meet technological requirements for hydrogen storage, owing to the weak nature of the H2 physisorption interaction with the substrate. Additionally, the maximum density of H2 that can be physisorbed is 70.58 kg m−3, the maximum hydrogen storage capacity of CNTs is thus limited by surface area constraints to 3.3 wt% (Züttel et al. 2004). Clearly layered materials must be modified if they are to be considered as realistic hydrogen storage materials, e.g. by external pressure (Gulseren et al. 2001) or doping (Shevlin & Guo 2006). It has recently been proposed that the addition of transition metal dopants can boost the hydrogen storage properties of carbon-based materials, by providing multiple sites for adsorption (boosting gravimetric storage), enhancing binding energy (allowing hydrogen storage at room temperature) and storing hydrogen in molecular form (increasing the speed of H2 loading and unloading). Electronic structure calculations on materials with Sc decoration on the fivefold faces of the fullerene C48B12 (Zhao et al. 2005), Ti decoration on one out of every two hexagons of an (8,0) nanotube (Yildirim & Ciraci 2005), and highly dispersed Ni atom decoration on a CNT (Lee et al. 2006) all show enhanced molecular H2 storage capacity (greater than or equal to 8 wt%). All of these systems are carbon-based, but is it possible that this approach may be transferable to other hexagonal net systems such as boron nitride, and indeed may be applicable to systems with ring structures other than pentagons and hexagons.
Before we discuss the details of the calculations, we shall first of all briefly discuss the mechanism by which H2 binds to the transition metal (TM) dopant. This is known as a Kubas interaction (Kubas 1988), which involves the donation of electronic density from the σ-state of the H2 molecule to an empty d-state of the TM atom, thereby leading to backdonation of the electronic density from a filled d-state of the TM atom into the σ*-state of the H2, causing a partial dissociation of the molecule and increasing the binding energy (figure 10). This is a three-centre, two-electron mechanism and is the σ-bond analogue of the Dewar interaction (Kubas 2001), which involves the donation of electronic density from the π-bond of a C−C bond to an empty d-state of the TM atom, and then backdonates electron density from an occupied d-state to the π*-state of the C−C bond (figure 10).
A good first step is to follow the rational design methodology of Zhao et al. (2005), who used cyclopentadiene rings as molecular substitutes for the pentagonal rings of the C60 molecule. If we are to adapt this to boron nitride, we need a molecule with a similar structure to the basic hexagonal net of h-BN. We therefore use, as a molecular representative of the basic building block of the h-BN net, borazine (B3N3H6), which is the boron nitride equivalent of benzene (C6H6). Indeed, as for benzene the borazine molecule possesses six π-electrons (Steiner et al. 2002) and thus we would expect a Dewar type coordination with TM atoms. The interaction of H2 with borazine is very similar to h-BN, i.e. in both cases simply a van der Waals interaction, and so we obtain similar binding energies of −0.02 eV/H2 for borazine and −0.04 eV/H2 for h-BN.
Structurally, we find that, of several different configurations, the lowest energy site for the TM atom is above the centre of the borazine molecule (figure 11a). The calculated binding energies for individual TM atoms (and magnetic moments in Bohr magnetons in brackets) are (in eV) from Sc to Fe: −2.50 (3.0), −1.56 (4.0), −1.64 (5.0), −0.02 (6.0), −0.12 (5.0) and −0.67 (4.0). As can be seen, the strongest binding is with scandium and titanium and the weakest with chromium and manganese. The poor bonding of Cr and Mn to borazine is because of the inability of the nitrogen atoms to donate to the low-valency TM atom coupled with the weak electron acceptance of the boron atoms (Bridgeman 1998). The TM binding energies are significantly weaker than for benzene where the carbon atoms may more easily accept electrons donated from the d-orbitals of the TM atom. Such weak binding means naked Cr and Mn atoms are not stable and will migrate across the molecule. TM atoms are less strongly bound to B3N3H6 than to C6H6, except the Sc dopant, which is more strongly bound to the borazine ring with the binding increased by 0.8 eV.
We then consider the addition of hydrogen, taking B3N3H6Sc as an example. The bonding of Sc to borazine is via a Dewar coordination involving the donation of a d-electron from the TM atom to borazine leaving two valence electrons to reside on the Sc atom. Hydrogen may bind to this Sc atom either by forming a dihydride with separate hydrogen atoms strongly bound to the Sc atom or as a dihydrogen ligand where the molecule binds directly with the Sc atom, as shown in figure 11b. We found that the dihydride was lower in energy by 1.04 eV over the dihydrogen ligand. The dihydride structure is preferred for all dopants, with energy differences between dihydride and dihydrogen becoming less from Sc to Fe, with the Cr and Fe dopants not possessing a dihydrogen structure at all. Addition of hydrogen does not change the bond lengths between borazine and the TM atom by more than approximately 0.1 Å. All the dihydride ligands (TM-2H) are bound to the borazine molecule, with binding energies of −1.00, −1.30, −1.05, −0.84, −0.30 and −1.58 eV/H2 for Sc, Ti, V, Cr, Mn and Fe, respectively (magnetic moments do not change). The addition of hydrogen stabilizes the Cr and Mn dopants on borazine owing to the hydrogen TM–H σ-bonds electronically tuning the TM atom to the ring, as examination of the electronic density of states shows that the valence states of the now hydrogenated TM atoms overlap the valence band of the borazine molecule. The hydrogenated Cr and Fe atoms now reside off the centre, forming strong bonds with two N atoms and one B atom.
In principle, the maximum number of hydrogen atoms that can be adsorbed on a TM atom should follow the 18-electron rule (Mingos 2004),where 18 is from the fully occupied s-, p-, and d-orbitals, nν is the valence electrons of the TM atom and 6 is the number of electrons contributed by borazine. Therefore, all the TM atoms should be able to adsorb additional H2 or 2H ligands (figure 11c). All successive hydrogen binding energies are shown in figure 12. Of note is that the second and subsequent H2 adsorptions are molecular in nature, meaning that the binding energy is substantially decreased.
As the overlap between the d-orbitals of the TM atom and the σ-bond of the hydrogen molecule increases with atomic mass, the hydrogen-binding energies tend to become larger. The Sc, Ti, and V atoms adsorb the largest number of H2, four H2 molecules each. Monohydride adsorption for large H2 concentrations was not observed, in contrast with the predictions of the 18-electron rule, which could be due to several reasons. The rule may be inappropriate for this system (Landis et al. 1998); or the weak aromaticity of borazine (about half that of benzene) may have an effect (Kiran et al. 2001). Of the three TM dopants, Sc binds too weakly to H2 and for low H2 capacity, V binds more strongly to H2 as well as being less stable on the borazine support. Therefore, the optimal TM dopant for this system is titanium.
Interestingly, boron nitride can form closed shell clusters that are composed of hexagonal and tetragonal BN rings (Shevlin et al. 2008) which is not the case for carbon fullerenes, which can only form low-energy structures composed of pentagonal and hexagonal carbon rings. The binding of TM atoms to tetragonal rings might lead to improved hydrogen storage characteristics, to test which we modelled 1,3-diazadiboretidine (B2N2H4), its interaction with a single Sc atom, and the hydrogenation properties of the composite system. The DFT calculated binding energy for the Sc atom to the ring is −2.31 eV, less than that for the borazine system. The Sc atom has a lower energy on inserting into the B−N bond, but the nature of the binding of the Sc atom to the B2N2H4 ring is the same as for the binding to the borazine ring (figure 13). We also find that upon hydrogenation, the composite system as a whole binds six H2 at ‘good’ binding energies (see start of section), with successive binding energies (for each H2) of: −1.37 eV/H2 (1H2), −0.26 eV/H2 (2H2), −0.53 eV/H2 (3H2), −0.41 eV/H2 (4H2), −0.23 eV/H2 (5H2) and −0.18 eV/H2(6H2). Clearly, TM doping of fourfold rings gives superior gravimetric storage performance than TM doping of sixfold rings.
Simulations of TM-doped carbon and boron nitride systems usually consider dopants that are dispersed uniformly and atomically on, for example, the surfaces of perfect CNTs. However, experiment shows that nickel forms approximately 1 nm diameter nanoparticles when adsorbed on CNTs (Lee et al. 2006), while simulation shows that for atomic Ti deposition on C60, the energetically lower structure involves nanoparticle coalescence rather than widely dispersed TM atom dopants (Sun et al. 2005). In order to minimize metal–metal interactions we need to ‘pin’ the dopants in place on the carbon nanotube. Vacancies in particular are of interest, as in the limit of low TM atom concentration (of the order a few per cent) the structures are intermediate between CNTs and metallocarbohedrenes (e.g. Ti14C13), of possible interest for hydrogen storage applications (Zhao et al. 2006).
As a case study, we consider the interaction of a Ti atom with a (8,0) carbon nanotube. As it has recently been shown by highly accurate quantum Monte Carlo simulations that the adsorption of hydrogen on CNTs (Okamoto & Miyamoto 2001) and light-element-doped fullerenes (Na-Phattalung et al. 2006) is more accurately described by the LDA than the GGA exchange correlation functional, the majority of our simulations were performed using the LDA, although several calculations were checked using both functionals. A dopant concentration of 0.78 per cent was modelled using an orthorhombic simulation cell of 20.8×20.8×6a Å3, where the last cell parameter is oriented along the tube axis and a is set to the C−C bond length of graphite, 1.41 Å for the LDA and 1.42 Å for the GGA.
A bare carbon vacancy on a carbon nanotube is unstable with respect to reconstruction to pentagon and nonagon rings (5,9), so we considered the Ti adsorption on this reconstructed defect structure. We examined several possible adsorption sites: (i) the centre of the fivefold ring, (ii) the centre of the ninefold ring; (iii) on the [5,6] junction, (iv) on the [5,9] junction, and (v) on the [9,6] junction. For all cases, the vacancy reconstructs such that the Ti atom inserts into the vacancy (figure 13). There are no physisorption states for adsorption on the vacancy defect itself. The adsorption of Ti on these reconstructed defects and insertion into the vacancy is barrierless in our simulations, indicating that in real systems the stabilized vacancy structures are easy to modify by Ti doping. The Ti resides above the nanotube by 1.15 Å. We calculated the binding energy of the Ti atom taking the reference configuration as the energy of a Ti atom in bulk hexagonal titanium (which is experimentally the more likely option for a reference state) and obtained a value of −1.61 eV. For this system coalescence of two Ti atoms in adjacent vacancies will not occur. A Bader partial charge analysis was performed to determine the magnitude of charge transfer from the Ti atom to the carbon support, where we find that the Ti atom donates 1.42e to the carbon nanotube. As opposed to Ti bound to a perfect CNT, the Ti dopant is in a different electrostatic environment. Therefore, the hydrogen adsorption properties of this substitutional dopant are expected to be different from those of an interstitial dopant, especially as it only binds to three carbon atoms.
We consider only exterior hydrogenation of the tube as adsorption occurs on the exterior in molecular form (Bauschlicher 2001). Additionally, if hydrogen is present in the interior of the tube it will chemisorb endohedrally to the carbon network and not be present in molecular form (Meregalli & Parrinello 2001). The geometries and the binding energies of the hydrogenated states are shown in figure 14. Initial hydrogen adsorption is on the Ti dopant and the carbon atom that forms a Ti−C bond directly along the tube axis (figure 14), as this results in the least stressed bond. This first hydrogenation is in atomic form; the dopant induces dissociation. We note that this does not cause the Ti−C bond to break, as the Ti dopant forms a σ-bond with the atom next to the carbon vacancy while the hydrogen atom binds to the π-band of the nanotube that is oriented perpendicular to the tube surface. Other structures, where one H is on the Ti and the other H is on the C2 or C3 atoms, are higher in energy by 0.44 eV. If the hydrogen atoms are placed on the exterior of the tube so as to attack the carbon atoms that bind the Ti, then the resulting structure is 0.13 eV higher in energy. It is thus unlikely that C−H bond formation will replace C−Ti bond formation. Initial hydrogenation of the Ti atom was found to be molecular in nature, but atomic adsorption on the Ti and C atoms was found to be lower in energy by 0.44 eV. However, subsequent hydrogenation is molecular in nature, thus implying good kinetics for hydrogen storage and release.
Both LDA and GGA calculations bind 5H2, with 4H2+H bound to the Ti and one H bound to the C atom. For the LDA, the H2 molecule bond length is stretched significantly (from 6 to 11%), whereas for the GGA, the H2 bond length is stretched only slightly (0–1%), which is reflected in the EBind; for the LDA simulation EBind varies from −0.34 to −0.63 eV/H2, all in the range of desired binding energies for technological applications. However for the GGA simulation this is not the case with EBind ranging from −0.05 to −0.46 eV/H2. The ‘true’ EBind is likely to be between the two values. We emphasize that, for both simulations, the binding of H2 to the Ti is enhanced with respect to the H2 binding to a pure carbon nanotube; doping with titanium enhances the thermodynamics of hydrogen storage. Finally, we note that the Ti atom directly binds 4H2+H, not the 5H2+H expected from the application of the 18-electron rule. We attribute this to the nature of the binding of the dopant to the carbon atoms, with the mixture of atomic orbitals involved preventing pure Kubas-type liganding. However the maximum number of H2 binding to Ti for this system can be formally described by a 16-electron rule, although the structures we obtain are not square-planar and so would seem to exclude this possibility.
Next, we consider the process of hydrogen storage and release. The change in free energy per metal atom upon successive hydrogenations is defined aswhere nH is the number of hydrogen atoms (and thus represents double the number of H2 molecules), nm is the number of metal atoms, and μH is the chemical potential of hydrogen, determined via the method of Reuter et al. (Reuter & Scheffler 2002) at T=298.15 K. The free energies as calculated with different exchange and correlation density functionals as a function of pressure are shown in figure 15. Full hydrogenation of the ‘high-c’ structure occurs at 58.6 atm (64.3 for the GGA). These pressures for hydrogenation are of the same order of magnitude as those required for the hydrogenation of nanostructured carbon materials, but with much better thermodynamic properties.
These calculations demonstrate that, by itself, boron nitride, regardless of the presence or absence of native point defects, cannot meet the current technological requirements for hydrogen storage applications. We have demonstrated that TM doping can improve both thermodynamic and gravimetric hydrogen storage properties as the chemistry of the B−N bond permits the formation of fourfold rings, BN clusters in principle are ideal for TM-assisted hydrogen storage. Finally, we have demonstrated that it is possible to create a nanostructured CNT so that it may pin Ti dopants, preventing coalescence and yet still binding a significant amount of H2 per dopant (5H2 per dopant, for a total gravimetric capacity of 7.1 wt%). Simulation may thus act as a way of predicting and designing novel hydrogen storage materials.
(c) Energy storage and conversion: from fuels to battery materials
Alternative means of storing, transporting and utilizing energy require a different set of functional materials as discussed in this section.
(i) Methanol synthesis via CO2 activation
As discussed above, one of the most efficient and environmentally friendly ways to store and use chemical energy is by using hydrogen, which has given rise to the prospect of the hydrogen economy. In the first part of this review, we have outlined the challenge to materials science presented by hydrogen production, focusing on water-splitting, but we note that hydrogen is also a by-product of a number of industrial chemical processes. The next challenge arises from the necessity to store and transport hydrogen to a site of consumption, and we have highlighted recent materials research in the field of hydrogen storage. Here, we consider another means of hydrogen storage and transport, methanol, which is already widely used as a fuel in the automobile industry and beyond, and is an alternative to more traditional fuels obtained by oil refinement.
For methanol or higher alcohols and related compounds, the challenge in computational materials science is not in the compound itself but rather in its production, or more particularly, the industrial catalyst. The first widely used catalyst was chromia (and alumina)-doped zinc oxide, until the introduction of the ICI low-temperature catalyst comprising copper supported on zinc oxide doped with alumina in 1965. The main advantage of the latter system is in the significantly lower pressure and temperature operating conditions required for efficient synthesis. With both catalysts, methanol is produced from feed gas (syngas), which is a mixture of H2, CO2 and CO. Thus, the production of methanol not only presents an attractive alternative to hydrogen storage, but also allows us to use carbon dioxide—one of the main anthropogenic factors in climate change. However, carbon dioxide is difficult to activate and our understanding of the catalytic processes involved at an atomistic level has been lacking until very recently.
In a series of studies (French et al. 2001, 2003, 2008; Bromley et al. 2003; Catlow et al. 2003, 2005, 2008; Sherwood et al. 2003), we have first concentrated on issues of catalysis on ZnO, which has been extensively studied experimentally, as similar reaction mechanisms were postulated for the two catalysts (Bowker et al. 1981; Bailey et al. 1995; Harikumar & Santra 1996; Nakamura et al. 1998), even though some details could be distinctly different. We have considered the adsorption of methanol precursors including carbon dioxide, formate and methoxy ions, and have calculated the equilibrium structure and binding energies of these and related species, especially relating to the following reactions:where the gas phase and adsorbate species are indicated by (g) and (a), respectively. By studying the energetics of these reactions on the oxygen-terminated surface of ZnO, we have confirmed surface interstitial vacant sites as the catalytically active sites on the surface of ZnO. The formation of oxygen and zinc interstitial surface vacant sites is inherent to the mechanism of stabilization of the two polar surfaces of ZnO, Zn-terminated (0001) and O-terminated (000); where about a quarter of terminating atoms have to be removed from the surface to cancel an unsustainable macroscopic dipole across the crystal. These sites present natural deep traps for electrons donated to the conduction band of ZnO by aluminium impurities as shown in figure 16. Carbon dioxide is then activated on adsorption when it acquires a negative charge and bends to be readily transformed into the key reaction intermediate, formate, on hydrogenation. Subsequent hydrogenation of the formate by pre- or co-adsorbed surface hydrogen leads to the eventual formation of methanol. On methanol desorption, an oxygen atom occupying the active site can then be removed either by reaction with CO or hydrogen, the latter may be either gas phase or pre-adsorbed on the surface. Curiously, as the surface oxygen is readily scavenged by carbon monoxide, the active sites could in fact be mobile on the surface, when the original active site could swap its role with a regular oxygen surface site that would become vacant as a result of this reaction.
The role of the active site and the main source of carbon in methanol synthesis over both the industrial catalyst Cu/ZnO/Al2O3 and the doped ZnO support have been often misunderstood. Experimental studies from the group of Waugh and others have shown unambiguously that CO2 is the key species that goes on to form methanol over both catalysts, while CO along with H2 forms and maintains active sites (Ueno et al. 1970, 1971; Bowker et al. 1981; Tabatabaei et al. 2006).
Having established the first model of the ZnO catalyst, we carried out our investigations into the formation of the metal-supported catalysts and have considered the nature of interaction between the metal and support as well as the modes of metal growth on the oxide support. From early on we have identified surface-vacant Zn sites (shown in figure 17) occupied by the metal such as copper or gold that serves as an anchor for further cluster growth, which allowed us to assess the energetics of the growth using binding and nucleation energies. This study has been particularly important as novel catalysts have been advanced in recent years that involve small sub-nanometric clusters of gold and other coinage metals (see Phala et al. 2005 and references therein). The new catalysts are much easier to control while their performance is comparable with that of the industrial catalyst. Here lies the hope of further optimization of the catalyst, perhaps leading to further reduction of operating temperature and pressure, longer life and higher turnover.
(ii) Modelling of solid oxide fuel cell electrolytes
Solid oxide fuel cells (SOFCs), illustrated schematically in figure 18, are playing an increasingly important role in energy-generation technologies. There are major issues relating to materials performance for both the electrodes and electrolytes. Here, we concentrate on the latter, highlighting the input that modelling has given to fundamental processes controlling ionic transport in these materials.
The basic requirement for an SOFC electrolyte is that it has high oxygen ion mobility; and despite many decades of study, the materials used in current SOFCs are fluorite-structured oxides—either yttrium-doped zirconia or gadolinium-doped ceria. (Note that ZrO2 is only stable in the fluorite structure with relatively high concentrations of low valence dopants.) The basic defect chemistry of these systems is simple. The low valence dopant (Y3+ or Gd3+) is compensated by the formation of oxygen vacancies; and in common with all fluorite-structured materials, anion vacancies are mobile, with both calculations and experiment giving a low value for the activation energy of approximately 0.4–0.5 eV in both zirconia and ceria (Butler et al. 1983; Balducci et al. 1997; Kamiya et al. 2001; Dong et al. 2004) for a simple migration process involving a linear path between two oxygen sites, as shown in figure 19. A crucial feature concerns the interaction between the dopant ion and the vacancy: as these have opposite effective charges, they will have an attractive interaction and can form complexes, the simplest of which is shown in figure 19. Early calculations using interatomic potentials (Butler et al. 1983) were able to estimate the strength of these interactions, which will reduce oxygen ion mobility. Interestingly, work on CeO2 showed that the strength of the interaction depends strongly on the ionic radius of the dopant, implying that the binding energy has an elastic, as well as an electrostatic, component. The minimum binding energy (of 0.26 eV), for the complex shown in figure 19, was obtained for Gd as dopant—a result that fully aligns with the use of Gd as the preferred dopant, when ceria is employed as the electrolyte. A number of studies have been reported on ZrO2, most recently by Xia et al. (2009). They also show that again the dopant acts as a trap for vacancies, but it appears from both calculations and experiment (Catlow et al. 1986) that the most stable site for the vacancy in Y/ZrO2 is commonly in the next nearest neighbour position with respect to the vacancy, rather than the nearest neighbour position shown in figure 19; although the calculations of Xia et al. suggest that the energy difference may be small and concentration dependent.
Xia et al. also examined the important question of dopant segregation to the surface, which is a common phenomenon in ceramic materials such as zirconia. They reported an interatomic potential-based modelling study of the surfaces of the material and established that, as is generally the case with fluorite-structured oxides, the (111) surface has the lowest energy, with the (110) surface being the next most stable. Their detailed examination of the energetics of yttrium–vacancy complexes at both surfaces established that there is a driving force (of approx. 0.5–0.7 eV) for dopant–vacancy segregation to the (111) surface, which the modelling studies predict will result in dopant enrichment of the surface to a depth of 5 Å. In contrast, no driving force is obtained for dopant segregation to the (110) surface. The segregation of dopants and vacancies to the surface of the material is expected to have significant effects on the performance of the material as a solid-state electrolyte.
Both experiment and modelling continue to explore alternative oxide materials, for example, the apatite-structured oxides discussed by Islam (2010) in this issue; and indeed neither ceria nor zirconia are ideal electrolytes as the oxygen mobility is only sufficiently high at elevated temperatures. Development of a material with higher O2− ion mobility allowing lower operating temperatures would be a major breakthrough for the field.
(iii) Lithium intercalation in cathode material
As discussed in the article of Islam (2010), solid-state batteries, shown in figure 20, are generally based on intercalation chemistry, usually involving metal oxides, involving the cell reaction2.11The classic system is LixCoO2, which is still the basis of many commercial solid-state battery systems. Several other oxides (including ternary systems) continue, however, to be investigated, since, as argued by Islam, there is a strong incentive to develop new materials including the olivine-structured LiFePO4 material. Here, we summarize work on two other classes of material: vanadium oxides and spinel-structured, mixed metal oxides, concentrating on the role of simulations in calculating cell voltages.
The first and possibly the most significant contribution that modelling methods can make to understanding intercalation reactions and their role in solid-state batteries is the calculation of the energy of the lithium intercalation reaction. If the anode process contributes little to the energy of the cell reaction, and entropy terms are small, these energies may then be used directly to calculate the cell voltage, E, via the Nernst relationship. Moreover, by calculating the energy as a function of x (in equation (2.11)), we can calculate E as a function of x. Standard DFT (with GGA functionals) has been widely used for this purpose. An illustrative example is given in table 8 for the case of V2O5, which was investigated by Braithwaite et al. (1999).
The calculated discharge curve is compared with experiment: reasonable agreement is obtained, although there is a systematic underestimation of the cell voltage by approximately 0.5 V which had been observed by other workers (Ceder et al. 1997). However, given the complexity of the system and the relatively low ‘resolution’ of the calculated variation of E with x, the agreement is gratifying. The calculations are also able to identify the intercalation sites within the channel in the V2O5 host material.
A more complex system, the spinel-structured LixCoyMn4−yO8, was also investigated by Braithwaite et al. (2000). Their work was stimulated by the observation that these materials, when used as cathodes, could lead to exceptionally high cell voltages of greater than 5 V. The result is of particular interest as neither binary manganese nor ternary oxides show such large voltages. Table 9 shows the results of calculated cell voltages for a variety of compositions. If we systematically correct the calculated value by the estimated one of 0.5 V above, we find good agreement with experiment. We note that the very high voltages observed experimentally are reproduced when Co is present in the mixed metal oxide. These voltages are associated with the reduction of Co; but the analysis of the variation in electronic structure with the degree of intercalation reveals a complex process with subtle variations in electron density on oxygen as well as cobalt and manganese.
A rather more exotic example is provided by the work of Kganyago et al. (2003) who modelled lithium intercalation in the Chevrel phase material MgxMo6S8. This detailed study was able to identify intercalation sites and again to predict cell voltages. Together with the other examples discussed here and the work summarized in the article of Islam (2010), it provides a good illustration of the predictive nature of simulations in solid-state electrochemistry.
(iv) Tuning the critical temperature of vanadium dioxide with tungsten
While V2O5 is the most studied vanadium oxide material for battery applications, closely related VO2, and indeed other oxides with vanadium in lower oxidation states, pure and doped, have been demonstrated as potential cathode materials. Moreover, vanadium dioxide exhibits intriguing physical behaviour and thus provides an exciting challenge to our fundamental understanding of the structure–property relationship in this class of material. In particular, vanadium dioxide provides a classic example of a metal–semiconductor transition. With a critical temperature (Tc) of 341 K (Rakotoniaina et al. 1993), which is tuneable over a wide range by doping with, for example, tungsten, pure and doped VO2 is potentially an ideal material for many different applications: thermochromics, temperature-sensing devices, modulators and polarizers of sub-millimetre-wavelength radiation, optical data storage media, and variable-reflectance mirrors (Yin et al. 1996). The transition temperature of VO2 decreases at a rate of 11 K per atom% of niobium or molybdenum and 26 K per atom% of tungsten (Tang et al. 1985). Thus, the transition temperature can be lowered to 0°C with, for example, 3 atom% of W (i.e. δ=0.03) in V1−δWδO2 (Sella et al. 1998). Regarding the application of infrared-active coatings, Guzman (2000) reported that the ability to tune Tc by doping makes tungsten-doped vanadium dioxide useful in energy-efficient or smart windows, which find special applications in the architectural, automotive and aerospace sectors. An energy saving surface coating composed of a vanadium dioxide sol-gel solution is particularly useful for the coating of large glass or plastic surfaces such as windows, walls and windscreens, including windshields. It allows sunlight (both visible and infrared) to penetrate through windows and heat the interior at lower temperatures. When the temperature rises, it switches (as a result of the phase change) to block infrared light from entering the interior while continuing to allow visible light to pass through, thereby reducing the heating of the interior.
For VO2, we have investigated both its phase transition (Woodley 2008) and, for the tungsten doped material, the local structure about tungsten (Netsianda et al. 2008). Both phases of VO2—a monoclinic and a tetragonal, or rutile phase—can be constructed from one-dimensional chains of edge-sharing VO6 octahedra. In the ‘ideal’ high temperature, rutile phase, there is one unique V–V distance, labelled DM in figure 21a, between two edge-sharing octahedra. In the low-temperature phase, there is a Peierls distortion along each chain so that two unique V–V distances, labelled DS and DL in figure 21b, alternate. In our chosen monoclinic unit cell, there are two distinct chains, one along an edge and another down the centre of the unit cell, and in each there is one short and one long V–V distance. Thus, figure 21b shows one of the four equivalent models of the monoclinic phase; switch the sequence of short–long in one, other or both to obtain the other equivalent structures. There is also a small displacement of V perpendicular to the chains, as movement of V along the chains will result in a shorter V–V distance between neighbouring chains. Using interatomic potentials, we have found two (degenerate) vibrational modes that become soft for the tetragonal phase; we choose a linear combination so that Peierls distortion occurs along half of the chains in the crystal (i.e. along one of the chains in our monoclinic cell) accompanied by smaller perpendicular distortions in the other chains, as shown in figure 21a. To investigate the nature of this phase transition, we modelled many Peierls distorted structures with the change in the fractional coordinates from the tetragonal phase defined by the magnitude of these two vibrational modes, (L1,L2). The lattice energy as a function of (L1,L2) is shown in figure 22a. There are four stable, or low-lying minima; four equivalent monoclinic phases, which are possible in our chosen monoclinic cell, (L,L), (−L,L), (L,−L) and (−L,−L). The free energy, which includes the zero point energy contribution from the other vibrational modes, is shown in figure 22b. There is a discontinuity in the free energy, which marks the boundary between the two phases and implies, by definition, that VO2 has a first-order phase transition. Thus, from this surface we can clearly identify an ergodic region associated with the high temperature phase centred at (0,0). Moreover, each of the low-temperature minima has split so that |L1|≠|L2|, i.e. the magnitude of the Peierls distortion for the low temperature phase is different for neighbouring chains. It would be interesting to repeat our simulations for a larger monoclinic cell, one which contains more than two chains, in order to see whether magnitude is different along each unique chain. It is important to note that there is a small barrier between pairs of minima, so both are within the same local ergodic region—so over an observation time, at low temperatures, we would observe |L1|=|L2|. As temperature increases, the decrease in the free energy surface, or ergodic region, is faster for the tetragonal phase than the monoclinic phase(s), and eventually, at approximately 300 K, the ergodic region at (0,0) contains the global minimum—see figures in Woodley (2008).
In order to model the tungsten-doped phases, we used a Mott–Littleton approach, where a quasi-continuum approximation is employed to ions that are outside a spherical region centred on the tungsten ion. Thus, we modelled an infinitely dilute concentration of dopants, i.e. did not consider interactions between dopants or defect complexes. The local structures about the tungsten ion are shown in figure 23, where we highlight a short cation–cation distance within a chain using a stick. Interestingly, for the high-temperature phase the chains relax to create a short V–V distance on either side of the tungsten, i.e. a localized Peierls distortion occurs in the chain containing the W ion. Also, there are no short W–V distances for the low-temperature phase, so the alternating sequence—short–long–short–long—is broken to create two defect centres, the tungsten ion and the vanadium ion that has two short V–V distances. From the calculated substitution energies, 2.344 and 1.963 eV for the low- and high-temperature phases respectively, we can estimate the change in the critical temperature upon doping,2.12where Tc is the critical temperature of the pure phase, dH is the change in enthalpy (internal or lattice energy) of the pure phases and n is the fraction of substituted sites; see Netsianda et al. (2008) for a derivation of equation (2.12). For a 1 per cent solution of W we predict a reduction of 48.4 K, whereas it is observed to be a reduction of 26 K (Tang et al. 1985).
In this section, we have discussed several materials for chemical energy storage and conversion. However, direct energy conversion from sunlight is another route of interest for clean energy storage and we discuss in the next section novel materials currently proposed for solar cells. A reverse process, in which electrical power is converted to light, is also of high interest and, perhaps, not surprisingly, solar cell materials play an important role here.
(d) Materials for light absorption and emission
We start our discussion with an exciting range of new ternary materials, then turn the focus on GaN, a traditional optoelectronic material that finds new applications in solar cells, and complete with a summary of defect properties of ZnO.
(i) Inorganic solar cells
The direct conversion of visible light to electricity is a process with the potential to offset our current dependence on fossil fuels. Solar cells facilitate this process, through the photovoltaic effect, by extracting an electrical voltage from a semiconducting material irradiated with light; the maximum voltage obtainable corresponds to the electronic band gap of the material. Forty-four years after Becquerel’s photoelectric experiments, Fritts built the first solar cell using Se in 1883, with a modest visible light to electricity conversion efficiency of 1 per cent (Fritts 1883). It was not until the 1950s that silicon solar cells were successfully fabricated and commercialized with 6 per cent efficiency. After over 50 years of research, efficiencies have risen to above 20 per cent and crystalline silicon has become synonymous with photovoltaic devices; however, silicon is expensive to process and owing to a low optical absorption coefficient, solar cells based on silicon are generally thick and expensive. Alternative semiconducting materials, which are strong absorbers of light (direct band gap), require less active material, resulting in cheaper ‘second-generation’ thin-film solar cell devices (Luque & Hegedus 2003). More recent activity into ‘third-generation’ solar cells has focused on going beyond the theoretical conversion limit of 30 per cent for a single junction device (Shah et al. 1995), and includes multi-junction devices containing more than one active absorber material, multiple-exciton generation from quantum-confined systems and intermediate band gap solar cells.
The two leading commercial thin-film solar cells, based on CdTe and Cu(In,Ga)Se2, respectively, can reach conversion efficiencies approaching 20 per cent. In both cases, the absorber material acts as the active p-type semiconductor, which is matched with a wider band gap n-type semiconductor such as CdS, to produce the nominal p–n heterojunction that enables charge-carrier separation; a schematic of a high-efficiency CdTe thin-film solar cell is shown in figure 24. To improve further the cost-to-efficiency ratio, one would like to replace the Se and Te anions and the Cd and In cations by cheaper and more sustainable (abundant) alternatives. In the following, we will discuss progress in our understanding of leading candidates towards sustainable solar cells suitable for widespread deployment.
Kesterite thin-film solar cell absorbers
Multi-ternary chalcogenide semiconductors can be formed by taking simple binary compounds and cross-substituting to form ternary and quaternary systems (Goodman 1958). This breeding process can be used to tune the structural and electronic properties, e.g. lattice constant and band gap, as with traditional III–V semiconductor alloys such as Ga1−xInxP. For example, taking the transition from ZnS to CuGaS2, where two Zn(II) species are replaced by Cu(I) and Ga(III), one step further by replacing two Ga(III) by Zn(II) and Sn(IV) results in Cu2ZnSnS4 (CZTS), a I2–II–IV–VI4 semiconductor. The same methodology can be employed to construct an entire class of stable quaternary semiconductor systems (Chen et al. 2009b), in which the tetrahedral networks satisfy the octet rule of Lewis (1916), as shown in figure 25.
CZTS has become the subject of intense interest because it is an ideal candidate absorber material for thin-film solar cells, which combines strong optical absorption with the optimal band gap of 1.5 eV for a single-junction device. In comparison to current thin-film solar cell materials, CZTS contains abundant elemental components, and it has the additional advantage of being adaptable to solid-state, electrochemical and chemical-solution based growth techniques, which can be used to tune the material’s morphology, stoichiometry and intrinsic p-type conductivity (Scragg et al. 2008). A record power conversion efficiency of 9.6 per cent was recently reported for CZTS solar cells (Todorov et al. in press). In previous work, we explored the competing ground-state phases of CZTS and demonstrated that alternate cation ordering only has a weak effect on the electronic properties (Chen et al. 2009a). The two competing ground-state structures, kesterite (space group ) and stannite (space group ), are based on a 1×1×2 tetragonal expansion of zinc-blende with different ordering in the II–IV cation sublattice. As with the chalcopyrite structure of CuInSe2, the kesterite mineral structure is favoured as it maximizes the Coulomb energy: based on the formal oxidation states and the ideal case of c/2a=1, the difference in electrostatic energy between the two structures is 32 meV atom−1.
One question arises as to the source of intrinsic defects in these I2–II–IV–VI4 materials and how this could be used to optimize the growth of more efficient CZTS solar cell devices. The Madelung site potentials in the kesterite and stannite structures are similar (table 10), so we expect that the cation ordering will not significantly affect the fundamental defect chemistry of these materials.
From DFT calculations on kesterite CZTS (Chen et al. 2010), we have identified two principal low-energy processes that can generate p-type carriers, which are related to copper deficiency and excess, respectively. This study was performed using the PW91 exchange-correlation functional with a 64-atom periodic supercell. In semiconductors containing Cu(I), the anticipated process is the formation of copper vacancies:where the ionization of the copper vacancy generates a hole in the valence band (),The competing process is the incorporation of excess copper on the Zn sites:where the ionization of this copper anti-site also generates one electron hole,Both point defects are easily formed and have shallow ionization levels, so that they may both actively contribute to the intrinsic electrical conductivity; however, analysis of the accessible chemical potential range, calculated from the relative stabilities of CZTS and the competing binary and ternary components, suggests that Cu-rich/Zn-poor conditions are required to avoid precipitation of ZnS and SnS (Chen et al. 2010). Under the equilibrium conditions required to form high-quality CZTS samples, we predict that CuZn will be the dominant acceptor species. If excess Zn is present, it will provide ionic compensation through formation of the inverse anti-site,which is an active electron donor,The acceptor and donor levels produced by the Cu and Zn anti-sites follow the character of the valence (anti-bonding Cu d−S p) and conduction (predominately Sn, Zn and S s) bands, respectively, as shown in figure 26.
One omnipresent issue in the related chalcopyrite (I–II–VI2) thin-film solar cell absorbers is the self-passivation of electrically active point defects through complexation, e.g. in CuInSe2 (Zhang et al. 1998). If produced in significant abundance, such complexes can deplete charge carriers, restricting device performance, or conversely, may form internal heterojuctions that aid charge-carrier separation. In the quaternary systems, owing to the increase in the number of elements and possible defects, compensating pairs are to be expected. In CZTS, pairs have a predicted binding energy of 2.22 eV, giving an overall formation energy of 0.21 eV. These pairs can be viewed as partial disorder in the Cu−Zn sublattice, which has been observed in neutron scattering experiments (Schorr et al. 2007). In comparison, disorder in the Zn−Sn (0.86 eV) or Cu−Sn (1.26 eV) sublattices is less favourable. The driving force for defect complexation comes from a combination of Coulombic attraction and charge compensation of oppositely charged defects. The only deep states introduced into the band gap are from intermixing between Cu and Sn owing to their large size and chemical (site potential) mismatch.
These calculations highlight the important role that theory can play in the characterization of new material systems, and have identified CuZn and VCu as two key intrinsic acceptor defects in Cu2ZnSnS4. Progress in utilizing CZTS and related kesterite materials in solar cells will require further experimental investigation into high-quality crystalline samples, with careful measurements of the defect signatures including carrier concentration as a function of growth conditions. Combined with modern simulation techniques, the optimization process towards more efficient low-cost solar cells could be rapidly realized.
(ii) New interatomic potential for modelling gallium nitride
In contrast to low-cost metal chalcogenide thin-film photovoltaics, solar cells based on III–V semiconductors serve a different purpose, and encompass high efficiency devices more suitable for space, military and concentrated solar applications. A triple-junction device containing layers of GaAs, Ga1−xInxP and InxGa1−xAs holds the current light to electricity conversion record of 40.8 per cent (Geisz et al. 2008). To date, the main use of III–V nitrides (AlN, GaN and InN) has been in light-emitting diodes, lasers and photodetectors, but their application to solar cells is underway (Jani et al. 2007; Zheng et al. 2009).
A solar cell based on III–V nitrides, consisting of GaN as both the p-type window layer (front contact) and n-type back contact, with In1−xGaxN as the active absorber layer, has demonstrated high cell voltages and conversion efficiencies (Jani et al. 2007; Zheng et al. 2009), and is shown schematically in figure 27. One of the main challenges is to avoid phase separation of the In1−xGaxN alloy; the addition of InN (Eg=0.7 eV) is required to reduce the band gap of GaN (Eg=3.5 eV) in order to harvest visible photons efficiently.
Addressing the mechanical issues associated with solar cell construction and degradation requires theoretical methods that allow for the simulation of large system sizes, for which classical interatomic potential functions are ideally suited. An appropriate pair potential model would provide a robust framework onto which further studies could be performed, including hybrid quantum mechanics/molecular mechanics simulations of defect and optical properties. Here, we present a new potential model for GaN, based on the formal ionic charges that can reproduce the structural, elastic and vibrational material properties of both the wurtzite and zinc-blende polytypes. The system is modelled using a multi-region pair-wise potential that includes polarization of nitrogen through the shell model (Dick & Overhauser 1958). In this approach, parameters are optimized separately for the interaction with nearest, next nearest neighbours and the remainder of the lattice. For each range of typical interatomic distances, a suitable simple analytical potential function is used, and the differently parametrized functions are tailored using buffer regions with monotonic polynomial functions. The Ga−N interaction is modelled with a Born–Mayer potential. To prevent the unphysical approach of Ga–Ga and Ga–N ions, a Lennard-Jones r−12-term is introduced between each pair. The N−N potential for the first neighbour shell is parametrized in the Morse form, while the interaction with the remainder of the lattice is modelled through a van der Waals dispersion r−6-term. All parameters are accessible online (http://www.dfrl.ucl.ac.uk/Potentials/).
GaN has three thermodynamically accessible phases: wurtzite (B4), zinc-blende (B3) and rock salt (B1), with the ground-state phase being the B4 structure. The total energy difference between the B3 and B4 polytypes is 17.8 meV atom−1, while the B1 phase is 573 meV higher in energy, as calculated in the GULP package using the interatomic potential detailed in table 11. The experimentally determined critical pressure for the transition from tetrahedral to octahedral coordination (B4 to B1 phase transition) varies between 37 and 52 GPa. Explicit calculation of the phase enthalpies as a function of pressure at T=0 K, using our potential model, predicts a transition at 54 GPa to the B1 phase in fair agreement.
The equilibrium (T=0 K, P=0 GPa) structural properties of each phase of GaN are listed in table 11. Good agreement is observed between the experimental properties and those calculated using our model potential, including the gamma point phonon frequencies. Table 11 also contains results from an earlier interatomic potential based on a partial charge model (Zapol et al. 1997); it can be seen that our potential overall improves on the results of the previous work, in particular for the dielectric properties an accurate description is provided, which is essential for a reliable description of the defect chemistry.
The four principal types of charged native defects were studied in the wurtzite phase of GaN using the Mott–Littleton embedded-cluster approach implemented in GULP (Gale & Rohl 2003): nitrogen and gallium interstitials and vacancies. While GaN can be extrinsically doped to give n-type or p-type behaviour (Tansley & Egan 1992), intrinsic defects may significantly modify the properties of the host material and their effects cannot be neglected. Defect formation energies for Schottky (isolated Ga and N vacancies) and Frenkel (isolated vacancy and interstitial) pairs are presented in table 12.
Owing to the closed packed structure, the cost of interstitial formation is high, resulting in the large formation energies for both anion and cation Frenkel pairs. The calculations suggest that Schottky formation will be the dominant form of intrinsic ionic disorder. Reasonably good agreement is found with the partial charge potential. The defect energies obtained here using interatomic potentials are consistently higher than those from electronic structure calculations using periodic boundary conditions (in different charge states), but all approaches yield Schottky disorder as the lowest energy process.
The charged nitrogen and gallium vacancies result in slight distortions of the lattice in the vicinity of the vacancy. Introducing extra nitrogen into the lattice leads to the formation of a split interstitial as shown in figure 28. At equilibrium, the gallium interstitial is positioned in the middle of the hexagonal channel without significant distortions of the neighbouring lattice.
The model we have developed provides a sound basis for future interatomic potential studies of the III–V nitrides (AlN, GaN, InN). In particular, it could enable explicit studies of the (GaN)1−x(InN)x solid-solution critical for solar cell and white-light light-emitting diode applications, including the thermodynamics of phase segregation and the defect chemistry of the disordered system. Potential models are also essential when modelling defects using embedded cluster (or QM/MM) techniques, the power of which has been demonstrated in several recent studies (Catlow et al. 2008).
(iii) Defect centres in zinc oxide
We complete our account of light absorbing and emitting materials by a brief account of our studies into defect formation and properties of zinc oxide, which is to date, perhaps, the most widely used functional oxide material. As described above, it is the core part of the industrial methanol catalyst, but no less importantly zinc oxide is used as a support and window layer in GaN solar cells, as a phosphor and, along with GaN, as a blue-light-emitting material in diodes and lasers. The key physical properties of zinc oxide are its piezoelectric behaviour, pronounced n-type conductivity of virtually all samples and strong near band-edge and broad-band luminescence covering the whole visible optical region. An efficient and reliable route to the synthesis of p-type zinc oxide is still lacking but is essential to realize the full potential of zinc oxide in optical and optoelectronic applications. As with all semiconducting materials, the major factors that control the optical and electrical behaviour of zinc oxide are point defects. We have recently explored the defect chemistry of another transparent metal oxide material In2O3 (Walsh et al. 2009b), and work is in progress on related compounds. Experimental evidence, however, that would characterize intrinsic defect species and identify them firmly with particular spectroscopic signatures is scarce. Hence, computational techniques make especially useful contributions to this field.
Here, we outline results from our hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster (Sherwood et al. 2003; Sokol et al. 2004) investigations into intrinsic defect centres and a number of important impurities (Sokol et al. 2007), which are complemented where appropriate by semiclassical Mott–Littleton calculations (MM using a set-up similar to that applied above to study defects in GaN). Our first set of results is the comparative energetics for the ionic disorder in this material summarized in table 13.
We calculate remarkably similar energies of formation using the two methods for the two most favourable defect formation processes: (i) the Schottky pair yielding charged zinc and oxygen vacancies in their formal oxidation states and (ii) the anion Frenkel pair, in which a neutral oxygen vacancy, or an F-centre, is counterbalanced by a neutral oxygen interstitial. The neutral oxygen interstitial takes the form of a split interstitial in a typical closed-shell singlet peroxide species. These calculations have been complemented by an MM investigation into the defect complexation (Catlow et al. 2008), where we found the lowest energy complex formed by the Schottky defect pair with the nearest neighbour zinc and oxygen vacant sites arranged along the c-axis of the wurtzite hexagonal lattice, which brings the energy of formation down to ca 1.90 eV atom−1. Even though this energy is relatively small, it is still too high to explain typical concentrations of charge carriers in zinc oxide, and we must conclude that as with other dense oxides, the surfaces, grain boundaries and linear defects are the main source of intrinsic point defects in this material.
Next, we present the diagram of defect donor levels in the band gap of zinc oxide, shown in figure 29. In this diagram, all energies have been calculated with respect to the bottom of the conduction band and show the energies required to ionize respective defects, or transfer an electron from a defect to the bottom of the conduction band. To obtain these energies we have performed many-electron density functional calculations and have used the total QM/MM energies for the many-electron systems in the respective oxidation states of interest. In particular, we have used the full Born–Haber cycle to obtain the vertical ionization energies and not relied on one-electron Kohn–Sham states from the component density functional calculations, which has been, in fact, common practice in periodic boundary calculations. This procedure is discussed in detail in Sokol et al. (2007). Moreover, the approach avoids entirely the problem of artificial interaction between charged defects, which is one of the largest sources of error in the latter approach.
— Under oxygen-poor conditions, the formation of oxygen vacancies could be compensated through the Schottky mechanism (zinc vacancy formation), which does not affect stoichiometry, but concurrently through the anion Frenkel pair formation, in which neutral oxygen atoms are stabilized as interstitials. The interstitial oxygen then could readily migrate to the surface and sublime, which gives a clear rationale for the commonly observed oxygen substoichiometry of ZnO samples.
— The neutral and singly positively charged Zn interstitial defect is responsible for E1 and E3 (majority) donor bands from electrical measurements.
— Zn vacancies are the majority acceptors, in agreement with experimental assignment (Tuomisto et al. 2003) based on positron annihilation spectroscopic and other studies. This defect was found to be stable in five charge states. Exciton recombination at this defect species in these states is a source of the main photoluminescence bands: ultraviolet (an acceptor level at 3.2 eV in a donor to acceptor photoluminescence, DAP, transition), green (a triplet level at 2.5 eV), and red (at 1.9–2.0 eV).
— The neutral O interstitial in a split-interstitial peroxide configuration (at 2.8 eV) should also contribute to blue and green luminescence by an exciton recombination and DAP transition from donor Zn interstitials.
— O vacancies cannot be a source of green luminescence, but could contribute to near-gap (UV) and red-orange luminescence bands (at 2.1 eV and below) via the exciton recombination mechanism.
— Cu substitutional impurity, which is stable in ZnO in two charge states, is an efficient electron scavenger. The singly negatively charged Cu impurity was proposed as an E4 donor. Calculated defect levels (at 2.7 and 0.55 eV) are in good agreement with experiment, which established Cu as a distinct source of green luminescence from ZnO.
— Nitrogen, despite persistent suggestions in the literature, is not a shallow acceptor in zinc oxide, but is found to lie about 0.9 eV above the valence band in its neutral charge state and thus be a possible source of the blue band exciton recombination. This defect would present an efficient trap for electron carriers (ca 1.45 eV deep).
Following this work providing the first comprehensive rationale for the main body of defect phenomena in zinc oxide, we have extended our studies to the other potential sources of conduction and optical activity in this material, which are complemented by the studies into the mechanisms of diffusion, results of which will be presented elsewhere.
(e) Exploiting the nanoscale for current and future materials
Physicochemical properties of interest depend not only on the nature of the compound but also on the atomic-level structure, which at the nanoscale can vary dramatically with size. As we will see in the case of nitride materials, the variations do not necessarily follow quantum-size effects according to the confinement (effective mass) model, as is commonly assumed; however, similar structures can be found for different systems in certain size ranges, allowing the prospect of tuning material performance to a required specification, for example, relating cluster size and morphology to the energies and intensities of photo-absorption and luminescence bands.
This section considers the structural properties of nanoscale materials, relevant to energy technologies, and concerns different cluster stoichiometries. We first discuss 1 : 1 compounds, i.e. zinc oxide and related materials, a range of 2 : 3 (group 13 sesquioxide) compounds, followed by 1 : 2 (group 4 dioxide) compounds. We complete our account with a report of an example multi-valence compound, ceria, in which different stoichiometries coexist.
(i) Clusters of compounds with 1 : 1 stoichiometry
Global optimization studies (Behrman et al. 1994; Woodley & Catlow 2008) have predicted the lowest energy structures for each cluster size for different semiconducting compounds of 1 : 1 stoichiometry, (MX)n, and found remarkable similarities. Different structural motifs can be found for the global minimum in four different size ranges. The smallest clusters are typically planar, whereas in the next range of cluster size, the dominant structural motif takes the form of spherical bubbles, where all atoms are three-coordinated. Then stable clusters are found where some atoms are four-coordinated, thus forming onion-like structures, before the structural motif starts to resemble cuts from bulk phases. One particularly stable cluster is found for n=12 (figure 30), a spheroid configuration, which is often called a sodalite cage as it is the basic tile of the sodalite structure (SOD). The relative stability of different (MX)n clusters can be characterized by the energy of formation per formula unit with respect to one unit (for a dimer, M−N). Considering the most stable structure for each n in the lower range, we previously found a local minimum at n=12 for ZnO, and the next local minimum, although much shallower, at n=16 (Al-Sunaidi et al. 2008). The sodalite cage at n=12 is the first appearance of a spherical (also referred to as fullerene) cage with high symmetry (Th): for zinc oxide, stability of small clusters is often associated with high symmetry. On addition or removal of more units, similar sized, n=11 and 13, stable spherical cages are not possible. However, larger sized stable cages that are similar to that of the SOD, i.e. with the Th point symmetry group and where the number of surface hexagons is maximized, are possible (Shevlin et al. 2008). The n=16 cluster is also a spherical cage, but has a different point symmetry group, Td. Similarly, there is a series of ever-larger cages that resemble the n=16 cluster and have Td, rather than Th, point symmetry (Shevlin et al. 2008). Both sets of cages can be viewed as being constructed from eight sheets, in which the ions are arranged in a honeycomb, or hexagonal, pattern (figure 30). The sheets are bowed to create octahedron-based cages where, as in the sodalite cage, the strain induced by the curvature of the cages is reduced by the inclusion of just six tetragons. Finally, there is the special case of the n=24 cluster, which possesses the O point symmetry group (figure 30).
Now, let us consider another class of compound, nitrides, which however exhibit many similar properties. In particular, bulk III–V nitrides are of interest as they are medium to wide band gap semiconductors that have great durability over a wide range of conditions. Thus, they are suitable for applications in medical treatments, information storage and more generally, optoelectronic devices. The opportunities offered by modification of bulk properties via nanostructuring are immense. One pertinent example is the expected, and in many cases observed, opening of the band gap with efficient stabilization of excitonic states leading to light emission and lasing. In particular, h-BN has recently been found to be a promising deep UV emitter (Watanabe et al. 2004), and lasing has been observed, for example, for GaN nanowires (Johnson et al. 2002). It has been shown that the binding of excitons, produced upon photoexcitation, is strongly dependent on dimensionality (Wirtz et al. 2006). Thus, it is of interest to determine the structure, thermodynamics and electronic properties of nitride clusters as a function of size and morphology. As BN is the best studied nitride, i.e. has the largest amount of experimental data, we will concentrate on this material but will provide results for the other nitrides as a comparison.
As discussed earlier, while the B−N bond is isoelectronic with the C−C bond, BN cages are not structurally similar to carbon fullerenes as, due to the heteropolar nature of the BN framework, it is not energetically favourable for B−B and N−N bonds to form. Whereas the C60 molecule forms 12 pentagons in curving the graphite sheet, the same structure in BN systems would create B−B or N−N bonds at the pentagonal site, rendering the fullerene unstable with respect to B2 or N2 loss. Simulation has found that the lowest energy cages need both six- and four-membered rings to form octahedral structures (Seifert et al. 1997), with the general rule that six tetragons are needed to close the cage (instead of 12 pentagons for the case of C60). Although the experimental production of BN clusters was first reported in 1994 (Banhart et al. 1994), the successful synthesis of single-layer or reduced-layer cages did not occur for several more years (Golberg et al. 1998). For clusters greater than 1 nm in diameter, there is experimental evidence for the octahedral cages (Golberg et al. 1998). Golberg et al. have also shown that the cages are formed from ‘reduced-layer’ fullerenes and that the smallest diameter shell (approx. 0.4 nm) is a (BN)12 cage, regarded as the smallest possible high-symmetry magic number cluster in Jensen (1993).
We have identified several different structural motifs for zinc oxide (and related II–VI material) nanoparticle structures including rings, double rings, tubes and spherical cages: each type being more stable for a particular range of n (Al-Sunaidi et al. 2008). As several of the III–V nitrides are known to have similar bulk structures, it is reasonable to consider that the low-energy stable configurations of the nitrides may have similar topology to that found for ZnO. Thus, by rescaling the interatomic distances by the ratio of the respective ionic radii, we obtained initial candidate structures for the investigation reported here as shown in figure 30. The optical properties of BN cage structures are expected to approach h-BN single sheet properties as the cages increase in size. In contrast, for small structures, the optical properties should differ considerably from the single sheet properties, for example, owing to confinement effects. To understand the dependence of properties on the size of clusters, it is therefore essential to establish the sheet properties. The recently revised experimental value of the band gap for h-BN was found to be 5.971 eV (Watanabe et al. 2004).
To determine the relative stability of clusters of different size, morphology and chemical composition, we compare their energies of formation per formula unit with respect to gas phase diatomic molecules, which will be referred to in the following discussion as the cohesion energy:
Using relevant thermodynamic data and a Born–Haber cycle, we were able to make comparisons with that measured. Note that, currently, there are no reliable experimental data for the gas-phase diatomic nitrides of Al, Ga and In. Therefore, for the sake of consistency, we made use of results reported in the literature based on high-level computational studies. The reported dissociation energies have been obtained using the standard thermal term, (3/2)RT. To obtain the experimental bulk values of the cohesive energy, we have used the standard enthalpies of formation at 298.15 K.
All calculated cohesive energies are plotted in figure 31. As expected, the cohesive energies increase in magnitude and eventually saturate with cluster size. For clarity in this graph, we excluded points n=3 and 6, thus concentrating on closed single-wall, tubular and octahedral cage clusters. The n=60 double-bubble configuration is a special case, as can be seen from the graph where its energies significantly deviate from the smooth variation with size. The interpolation lines are in fact parabolic, which allows us to extrapolate the cohesion energies for n to infinity:We find remarkably good agreement of results obtained by two different simulation codes, VASP and DMol3.
It is clear that the strong linear terms determine the slope in these graphs for larger values of n, and we need to consider how this behaviour is related to the cage geometry. The cage size is characterized by a diameter, d, defined as the distance between the two most remote atoms in the cluster. To clarify further the dependence of the cage diameter on n and allow for extrapolation to , we plot in figure 32 the inverse square diameter as a function of inverse size, starting with n=36 as smaller clusters do not yield new insights.
For all four nitrides, the points fit very well on the interpolation lines (with R2>0.999), which have resulted from fitting to a parabolic function:The s/t ratio in this formula determines the value of n at which the quadratic term takes over from the linear, which proves to be about 0.09 for all four compounds. As seen in figure 32, the slope of the graph is dominated by the linear term at . We conclude that the asymptotic behaviour of the cohesion energy is described bywhich, in this form, coincides with the induced curvature energy, or bending energy owing to the strain within the surface of the cluster. Curiously, for spherical shell clusters, the total curvature energy depends on the compound and not on cluster diameter or size. Alternatively, we can consider that the total cluster energy could be split into the contribution of an idealized spherical shell (which, of course, could be realized only in the limit of infinite n) and the contribution from the energy of formation for six non-interacting tetragonal defects in the hexagonal sheet. Both possibilities are under investigation.
Finally, we turn our attention to the calculated one-electron energies of the nitride clusters, which provides our first insight into their optical properties. We have also plotted the HOMO–LUMO energies for the cage clusters as a function of 1/n in figure 33. We note that, starting with the n=36 octahedral clusters, the band gap values behave monotonically, with data lying tightly on parabolic interpolation lines with dominant constant and linear terms. The interpolation lines for AlN and GaN show the steepest negative slope, the BN line has a moderate positive slope while the InN energies are practically flat throughout the series of single-wall clusters. For BN, the band gap value extrapolated to is very close to the infinite sheet value obtained using the VASP code (4.69 versus 4.75 eV). Thus, we can conclude that the calculated trends do not concur with the currently widely held belief in the quantum confinement theory, which predicts strong size dependency of the optical properties irrespective of the actual nature of the material. Boron nitride displays the strongest variation in structure–property relationships for the range 12≤n≤36, whereas GaN and AlN only display a weakly monotonic behaviour of band gap for this size range, and InN has very similar gaps for all cluster sizes. Thus, nanostructuring in this range of cluster sizes is expected to be most effective for BN.
Silicon carbide clusters
The line of reasoning pursued in the studies outlined above could readily be extended to other compounds of 1 : 1 stoichiometry. Here, we will consider one such example, silicon carbide, a IV–IV compound, which is a prospective candidate for applications in energy, electronics and medicine. Of particular interest are microporous (Ciora et al. 2004) and mesoporous (Yang et al. 2006) materials, which exhibit chemical stability and mechanical strength along with relatively narrow pore size distributions. One route to the production of materials with porous architectures is by imitation of naturally occurring minerals with framework structure using differing chemical elements or whole molecular units as building blocks (Yaghi et al. 2000; Woodley & Catlow 2008), e.g. zeolites, as was recently proposed for nanoporous MgO and ZnO. In order to understand the properties of these framework units, it is essential to understand the basic properties, e.g. optical properties, of the cluster systems. Previous computational studies of SiC nanoparticles have usually been based on clusters carved from known dense bulk phases of SiC (Peng et al. 2007), and are of diameters that are similar to those that have been typically prepared experimentally, ca 4 nm. Small clusters, which consist only of a few atoms, were also investigated, with the focus on structure and properties of non-stoichiometric particles formed in early stages of gas-phase condensation (Erhart & Albe 2005). The latter has often been extrapolated to larger, medium-sized particles, where fullerene-like structures have been reported, based on pure and silicon-doped carbon or carbon-doped silicon clusters (Huda & Ray 2004).
Physically and structurally, SiC shows behaviour similar to BN as it forms a large number of low-energy polytypes based on the stacking sequences of hexagonal (honeycomb) sheets. The structural parameters of SiC are, however, closer to those of ZnO, which hypothetically could also exhibit polytypism even though only two such polymorphs have been observed so far, wurtzite and zinc-blende. ZnO is usually thought of as a more ionic material than BN, which favours stronger links between the hexagonal sheets, although the ratio of ionic radii between cations and anions in the material must also play an important role. The propensity of SiC to form the wurtzite structure (α-SiC), while the hexagonal BN-like structure does not seem to be stable, places it in this respect in between BN and ZnO, with a large body of experimental and theoretical work focused on the polytypism in this material (see Bechstedt et al. 1997). Therefore, it seems plausible that nanoscale cluster structures would be similar to those of the III–V and II–VI compounds. Indeed, we found from our calculations that even though clusters with point symmetry groups Th and Td were global minima for ionic systems, these structures are also highly transferable to other main-group compounds (see above). This is found also to be true for SiC, where these clusters are also global minima when compared with other potential cluster structures (Watkins et al. 2009). This is not a great surprise, as it has been shown that zinc-blende silicon carbide itself possesses substantial ionic character (Sabisch et al. 1995). The optical properties of these clusters were calculated so as to provide data that may be measured experimentally.
The optical properties of these clusters were calculated within time-dependent density functional theory (TDDFT) using the NWChem code, with the B3LYP exchange correlation functional modified using the asymptotically corrected functional of Casida and Salahub (Casida et al. 1998). A double zeta valence basis set was used. The gap for all clusters considered is given in table 14. With the exception of the (SiC)16 cluster, the gaps increase with increasing cluster size.
The excitation spectra of all of these clusters are found to possess many near-degeneracies. For example, the first six excitation states for (SiC)12 are found to lie within an energy range of 0.1 eV, as would be expected owing to the high structural symmetry. Therefore, we calculate and plot the density of states and predicted adsorption spectra up to 5.40 eV for the (SiC)12 and (SiC)16 clusters and up to 4.8 eV for the (SiC)24 cluster, as shown in figure 34. This procedure involves the calculation of 129 bound states for the (SiC)12 cluster and 200 bound states for the (SiC)16 and (SiC)24 clusters. Although we cannot compare our calculated spectra directly with experiment, we can compare with optical properties of nanocrystalline silicon carbide films grown by chemical vapour deposition (Rajagopalan et al. 2003), where it was found that high-quality films composed of nanocrystallites of size 2–10 nm with minimal hydrogen contamination have optical gaps of approximately 4.5 eV (Rajagopalan et al. 2003). Although the (SiC)12, (SiC)16 and (SiC)24 clusters are smaller than the observed nanocrystallite diameters, and although TDDFT calculation of excited states is known to underestimate energies by approximately 0.5 eV, the first states with significant oscillator strength occur in the range of 3.8–4.5 eV. Our calculated spectra for these clusters thus accord with available experimental data.
(ii) Plausible configurations of small sesquioxide (M2O3) nanoparticles
In a continuation of our search for low-energy small clusters of binary insulating and semiconducting compounds (Hamad et al. 2005; Al-Sunaidi et al. 2008; Woodley & Catlow 2008, 2009), we have investigated and report here a range of stoichiometric sesquioxide nanoparticles of M2O3 composition. Again, a hybrid genetic algorithm (Hamad et al. 2005) has been used to predict the five most stable configurations of (M2O3)n, for each value of M and n, where M=B, Al, Ga, In, Tl and Ce, and n=1–8. The materials concerned span a wide range of physical and chemical properties and are utilized as catalysts, lasing and more generally light-emitting materials, elements of complex photoelectric devices, etc. By exploiting the size–property relationship, we can further enhance these useful properties.
As the first step, we concentrated on the structure prediction for the smallest nanosized particles of sesquioxides. During the search a cost function that measures how good a structure is relative to others is required; we typically choose the energy of formation obtained using interatomic potentials. A range of interatomic potential parameters, available online (http://www.dfrl.ucl.ac.uk/Potentials/), were initially tried for the smallest clusters, n=1. Then with appropriate modification, just one set of potential parameters was employed for larger clusters. Here, we have employed a rigid ion model based on formal charges, i.e. +3 for the cations and −2 for the anions. The standard short-range interactions have been complemented by a strong repulsive term B/r12, where r is an interatomic distance and B is large enough to separate overlapping ions while having negligible effect for sensible r.
For (Al2O3)1, we have tried three sets of potential parameters (Bush et al. 1994; Minervini et al. 1999; Pandey et al. 1999). Similarly, for (Ga2O3)1, we tried parameters from a range of available potentials (Bush et al. 1994; Khan et al. 1998; Minervini et al. 1999). Finally, parameters for (In2O3)1 and (Ce2O3)1 were taken from the work of Minervini et al. With the majority of simulations, we predict the atomic ground-state structure to be a linear chain, shown in figure 35 along with the next two most stable configurations. Three exceptions, however, should be noted. (Ce2O3)1 adopts the dense configuration (figure 35c) as its ground state, which is also found for (Al2O3)1 using Pandey’s parameters and for (Ga2O3)1 with Khan’s parameters. The next most stable configuration for (Ce2O3)1 and (Al2O3)1 (using Pandey’s parameters) is linear. Even if the energy landscapes calculated with different potentials are different, these runs are still essential in order to generate a small set of plausible configurations. Therefore, we filter out plausible structures from the infinite number of possible, and so reduce the search space probed using ab initio approaches during further refinement.
The order of stability between local minima is expected to depend on the cation size. To study this effect, we have chosen one parameter set (Minervini et al. 1999), which is most complete and lacks only B and Tl; the missing data, however, were obtained readily by extrapolation. The full set of potential parameters is given elsewhere (Woodley 2007).
The difference in the formation energies of the linear and dense clusters for n=1, as well as the predicted bond lengths and angles, is shown in figure 36. While all interatomic distances increase with cation size, as should be expected, the bond lengths to single coordinated anions are much shorter than to doubly coordinated anions; and the dense configuration has the longest bond lengths. From the calculated energies, compounds with the larger cations are more stable in dense configurations. A linear extrapolation from Al to Tl would predict a structural transition point at the ionic radius of 108 pm; however, the change in the bond angle and differences in energy show a more complex picture.
Turning to larger clusters, the same n=2 ground-state configuration is predicted for all compounds. As part of the global search, aimed at locating the ground state, many local minima on the cost function hypersurface are found, and the corresponding metastable structures, for the lower energy configurations, are shown in figure 37. Our model predicts that, in general, the relative stabilities of each configuration do not change with cation size. The differences in stability between the predicted ground state (a) and the other investigated configurations (c)–(g) increase when the cation size is reduced. The second most stable configuration (b) of (Ce2O3)2 has a lower energy than that expected (if extrapolated from figure 37) for a group 13 cation of a similar ionic radius. For the smallest cation, boron, configuration (d) is less stable than might be expected if extrapolated. For larger cations, configuration (b) is very close in energy to (a), while for smaller cations configuration (b) relaxes to give configuration (a).
Next, we focus on the global minimum configurations—choosing alumina for a case study—as shown in figure 38. (Al2O3)3 resembles a pointed hat (shown in figure 38c with the anion point directed into the page), whereas (Al2O3)4 is a bubble and (Al2O3)8 adopts a tetrahedral shape. However, the sides, or faces, of the (Al2O3)8 tetrahedron are blown outwards as it is filled with two oxygen anions, each placed on one of six equivalent sites. The global minimum configuration of (Al2O3)2 forms a natural building block, or motif, that can be seen in the larger structures, particularly for n=4, 7 and 8.
Considering other compounds, five of the eight global minimum configurations of (B2O3)n, n=3–7 (figure 39a–e), are predicted to be different from alumina. The n=3, 4 and 5 configurations are the first three possible regular right prisms; two regular polygon faces connected by n squares, which is simply a cube for n=4. As shown in figure 39c, the cations form the corners of the prism (or corners of the polygon faces) and the anions are at the midpoints along the edges (or sides of the polygon faces) that are slightly bowed outwards. The next regular right prism is 1.76 eV higher in energy than the global minimum structure, indicating the relative instability of the six-membered rings formed by boron oxide.
As the ionic radius of the cation increases, from aluminium to thallium, the predicted global minimum configurations remain the same, with the exception of the two less symmetrical structures, n=5 and 6. The (Ga2O3)5, (Ga2O3)6, (In2O3)5, (In2O3)6, (Tl2O3)5 and (Tl2O3)6 structures are shown in figure 39 as configurations k, m, l, m, l and m, respectively. Furthermore, five of the eight predicted configurations for ceria differ from alumina, as seen in figure 39f–j. The insensitivity of the n=2, 4 and 8 configurations to the cationic radius and interatomic potentials used implies that these configurations are more likely to be predicted using an ab initio approach or, perhaps, adopted in nature. These configurations, however, could still change for smaller gap materials, where electronic terms were not properly accounted for by the current atomistic model.
(iii) Titania, zirconia and hafnia nanoclusters (MO2)n
So far, we have first discussed small clusters formed from 1 : 1 compounds, where bubble structures are typically the low-energy structures. It should be noted that there are 1 : 1 compounds that prefer to adopt structures that resemble cuts from the bulk phase, rock salt. The energy of formation is much lower for stoichiometric cuts that are perfect cuboids. Interestingly, the sizes, ncuboid, where perfect cuboids are possible matches the sizes, nmagic, that are more abundant in a plume of small particles created by laser ablation of a surface (Whetten 1993). Here, after addressing 2 : 3 compounds, we now turn our attention to a different stoichiometry, 1 : 2 compounds MO2, where M is a group 4 transition metal. Titania materials have been studied extensively owing to their unique photocatalytic properties (Thompson & Yates 2006; Chen et al. 2007) and other applications (Asahi et al. 2001; Gratzel 2001; Khan et al. 2002). The work on nanostructured TiO2 in particular is gaining importance (Kwon et al. 2008), and computational studies on TiO2 clusters provide a valuable insight into its structure and properties. Zirconia also has important applications because of its high melting point, low thermal conduction and high ionic conductivity. Its durability permits its use as a container for the safe disposal of high-level nuclear waste (Lumpkin 1999). Hafnia is increasingly studied, although it has the least current applications. Traditionally, only crystals of the order of a micrometre or larger were used, but in the last few years, there has been a considerable interest in the study of MO2 nanocrystals (Gratzel 2001; Frank et al. 2004; Kopidakis et al. 2005; Manera et al. 2005; Qiao & McLeskey 2005; Reddy & Khan 2005). One of the most promising applications of zirconia and hafnia is in the field of nanoelectronics. Owing to the increased miniaturization of transistors, the silica gates are reaching the lower size limit at which quantum tunnelling effects will cause failure, and therefore there is a strong incentive to find substitutes for SiO2. The high dielectric constant of ZrO2 and HfO2 makes them the most likely replacement (Wilk et al. 2001; Fiorentini & Gulleri 2002; Pasquarello & Stoneham 2005; Stoneham et al. 2005). We now review our work on the three compounds. A more detailed review of the computational studies of TiO2 clusters can be found elsewhere (Catlow et al. 2010). Although we use a similar approach to that discussed above for obtaining low-energy structures, one key point is that the influence of the d-electrons with regard to bond angles is not usually accounted for in our models based on interatomic potentials (IPs), i.e. when searching for plausible structures. Typically, this influence has a greater effect on the bond angles about under-coordinated atoms on the surface, which can always be corrected for in a later stage, e.g. further refinements of the plausible structures using DFT approaches.
In our initial study (Hamad et al. 2005), we applied several global minimization techniques to generate plausible structures of (TiO2)n clusters (n=1–15). In particular, we have employed a combination of Monte Carlo basin hopping and simulated annealing molecular dynamics, as well as different evolutionary algorithms, to search the landscapes defined using IPs. The plausible, or low-energy, structures were subsequently refined using DFT, and are shown in figure 40. On refinement, the change in the average Ti−O bond lengths decreases with cluster size from 0.074 to 0.030 Å, which is attributed to the fact that the IPs are fitted to reproduce the bulk properties of TiO2. The average Ti−O distance increases with cluster size, from 1.728 to 1.946 Å. The small clusters are isostructural to those predicted for silica (Catlow et al. 2010) and characterized by Ti ions concentrating in the core region of the cluster decorated by dangling, or singly coordinated, oxygen ions at the periphery. However, when the energy of formation is calculated using a DFT approach, after relaxation, the dangling oxygen atoms for titania do not remain in the plane containing the nearest Ti2O2 tetragon. For the larger clusters, the structures are compact, with a central octahedral and a surrounding layer of fourfold- and fivefold-coordinated Ti atoms, although there seems to be some energy penalty for particles containing the fivefold-coordinated metal atoms with a square base pyramid structure and dangling oxygen atoms. Here, a conventional double-bond notation is sometimes used between the dangling oxygen and the titanium, but this species is more likely to be a charge-transfer complex similar to that revealed recently for the silanone species in silica clusters (Zwijnenburg et al. 2009). The IP calculations tend to underestimate the dielectric screening of the charges localized on the dangling bonds and, thus, result in some preference for the more compact structures.
As titanium, zirconium and hafnium are from the same group of the Periodic Table and Ti(IV), Zr(IV) and Hf(IV) are chemically similar and isovalent, it is assumed that their respective configurational spaces are similar, and therefore the relative stabilities and structures of the more stable configurations will be similar. Moreover, isostructural configurations with the same order of stability for clusters of the oxides were found in an exhaustive search of stable configurations for n<4. Thus, in our study of zirconia, we started our main set of DFT optimizations for (ZrO2)n using scaled versions of (TiO2)n structures, and the similar procedure was followed for (HfO2)n (Woodley et al. 2006). Titania clusters tended to relax so that titanium atoms were in a more bulk-like environment, thus resulting in the creation of dangling oxygen atoms and the relatively constant value for ELUMO (discussed below) for clusters greater than n=3. Two of the minimum energy configurations for zirconia and hafnia, namely n=5 and 7, are more stable with a more densely packed configuration, which contained fewer or no dangling oxygen atom(s). There are similar trends, with respect to cluster size, in the electronic properties for the three compounds, especially between zirconia and hafnia. The predicted energies of the highest occupied molecular orbitals, EHOMO, lowest unoccupied molecular orbitals, ELUMO, and their differences for our lowest energy clusters are shown in figure 41. As could be expected for all ionic oxides, typically the HOMO is distributed mainly about the lowest coordinated oxygen anions and the LUMO about the highest coordinated cations. Thus, the HOMO is easily located when the configuration contains dangling oxygen atoms. One exception to this rule is found for (MO2)8, where the HOMO is distributed (as a direct result of the Madelung potential) over the lone pairs of the three oxygen atoms with coordination numbers 2, 3 and 3, rather than across all the two-coordinated anions. The variation in ELUMO, as a function of n, for (HfO2)n is mirrored, although always slightly lower in energy, by that predicted for (ZrO2)n. As the structures of (HfO2)n are similar, both in shape and size, to their (ZrO2)n counterparts, it is not surprising that ELUMO and the predicted energies for the highest occupied molecular orbital, EHOMO, and therefore their differences, are also similar. For the titania clusters, it is perhaps the very different variation in ELUMO that is more interesting; it is quite constant above n=3 as there is less change in the environment about the highest coordinated Ti, compared with Zr and Hf. The observed band gap for bulk titania and bulk hafnia is 3.75 eV (direct; Cronemeyer 1952) and 5.68 eV (Balog et al. 1977), respectively, whereas for bulk zirconia it is generally quoted to be about 5 eV (Navio et al. 2001). From quantum-size effects, one may expect that ELUMO–EHOMO would decrease as the cluster size grows, and perhaps this effect can already be seen for titania (an oscillation superimposed on a decrease from n=4). However, it is clear from figure 41 that ELUMO and EHOMO for zirconia and hafnia diverge with size, showing that the difference in the geometries of the growing cluster has more of an influence than the relative size, and the amplitude of oscillation found for titania is still significant. Moreover, although bulk-like features may be identified in the structure of our larger clusters, the fact that the band gap is largest for bulk hafnia whereas our clusters always have ELUMO–EHOMO larger for zirconia, it is clear that the electronic properties of such small clusters do not yet resemble those of the bulk.
We also computed the harmonic vibrational frequencies and infrared spectra for all three compounds (Hamad et al. 2005); two example infrared spectra are shown in figure 42. Computing infrared spectra, based on the DFT-optimized clusters, proved expensive (Hamad et al. 2005), thus we only considered n=1–8. From the computed harmonic frequencies of vibration, it is clear that some of the titania clusters have one or two frequencies with a much higher value (above 990 cm−1). These modes correspond to the vibration of the dangling oxygen atoms: all have two high frequencies corresponding to a symmetrical and asymmetrical stretching oscillation of the two dangling oxygen atoms except n=6 and 8, which have one and zero high-frequency/dangling oxygen atoms, respectively. The connection between high-frequency modes and dangling oxygen atoms can also be made for zirconia and hafnia. The stretching or compression of the bond(s), which connects dangling oxygen atom(s) to the cluster, results in a change in the cluster’s dipole moment and therefore an infrared active mode. Although the intensity of the infrared peak for dangling oxygen atoms will in general be smaller than that for the ‘bulk-like’ oxygen atoms oscillating in phase, this peak will be at a higher frequency as, for example, seen in only one of the infrared spectra shown in figure 42. The additional higher frequency peak provides a way of identifying the existence of dangling oxygen atoms.
(iv) Nanoclusters with multiple oxidation states: example of ceria
We conclude by considering ceria, which unlike zinc can readily adopt more than one oxidation state. Thus, for neutrally charged clusters, there is a greater range of stoichiometry to consider when the oxidation state of the cation is not fixed. Indeed, it has been demonstrated experimentally that below a certain size (e.g. 1.5 nm in the reported samples), ceria rapidly loses oxygen and gradually transforms from a dioxide to a sesquioxide form. This phenomenon manifests itself through anomalous lattice expansion (where average bond length within particles increases significantly; Tsunekawa et al. 1999, 2000; Wu et al. 2004). As discussed above, ceria is widely used in heterogeneous catalysis owing to its ability to uptake and release oxygen while preserving its crystal structure; most importantly, it catalyses oxidation of carbon monoxide to dioxide. The oxygen storage capacity of ceria strongly depends on the size and shape of the nanoparticles. At the lower end of the nanosize range, there have been only a limited number of investigations (Cordatos et al. 1996; Aubriet et al. 2009; Bromley et al. 2009).
In our investigation, as with the previous compounds, we have employed an evolutionary algorithm to search a range of landscapes, this time CexOy, as defined using a rigid ion model (two-body interatomic potential parameters are taken from a parametrization of a shell model used to simulate the bulk phase (Sayle et al. 1995), but using formal charges, i.e. +4, +3 and −2 point charges on Ce(IV), Ce(III) and O). For each landscape, we searched a number of times. Each was initialized with a different population of 30 random structures; then, for 200 cycles, a further 10 random structures and 30 offsprings were generated that competed with the current population to survive and become of the next population of 30 non-identical structures. Throughout the search, all new configurations are contained within a sphere with a radius of 8 Å and immediately relaxed so that during competition, comparison is made between energies of local minima. All structures in the last population are further refined using the shell model, and the lowest energy configurations found for small neutral clusters of CexOy (y<13) are shown in figure 43.
For clusters containing only cerium ions that have an oxidation state of IV, i.e. (CeO2)n, we show six structures. The configurations for the smallest three are isostructural to those predicted to be the global minima for silica (Catlow et al. 2010). However, when the missing electronic effects are included, the clusters can further relax, as reported earlier for titania, so that the bond angle for the n=1 cluster is less than 180° and likewise the singly coordinated oxygen atoms move out of the plane containing the nearest tetragon for n=2 and 3. DFT calculations reported by Chen et al. (2008) predict that the global minima for n=1–4 are all isostructural to titania; moreover, our lowest energy structures are also for sizes n=5 and 6. For clusters containing only cerium ions that have an oxidation state of III, i.e. (Ce2O3)n, we show four structures. The configurations for the smallest two agree with our earlier predictions; the hat and bubble structures predicted earlier for n=3 and 4 are not as stable after the configurations in the final population are refined using the shell model. In fact, we obtained a very symmetric cluster (Oh point symmetry)—a Ce6O8 octahedral bubble-like cluster with an additional oxygen atom within the centre. The Ce6O8 octahedral cluster is the smallest of the octahedral clusters that expose O-terminated facets—see Loschen et al. (2008) who also modelled, using a DFT approach, cation ordering (Migani et al. 2009) and oxygen vacancies (Loschen et al. 2008) in larger examples, Ce19O32, Ce44O80 and Ce85O160. The remaining clusters shown in figure 43 have mixed III/IV oxidation states, where typically the Ce(IV) is predicted to have a higher coordination than Ce(III). We repeated the search, for these stoichiometries, under the assumption that the additional electron on the Ce(III) ions is delocalized or shared between the cerium ions; an average point charge of 2y/x |e| for the cerium ions and Ce−O potential parameters were used. The lowest energy structure found for each size is shown in figure 44.
The lowest energy structures for Ce5O9, Ce6O11 and Ce7O11 are essentially the same; applying small atomic displacements, the more symmetrical clusters in figure 44 can be mapped onto those in figure 43 (note that the view of clusters, particularly for Ce5O9 and Ce7O11, in the two figures are different). For each of the remaining clusters, Ce3O5, Ce4O7, Ce5O8, Ce5O10 and Ce7O12, we predict a different lowest energy configuration, and typically the structures have higher symmetry when delocalized electrons are modelled.
In this survey of the nanoclusters, we emphasized the task of obtaining the atomic structure adopted by the materials of interest, but, as shown in our investigation into 1 : 1 bubbles, nanostructuring broadens the field of possibilities for modifying different physical and chemical properties of these materials for the use in energy technologies.
In this review, we have attempted to illustrate the wide ranging applications and impact of computational techniques in the field of energy materials. It is clear that these methods are now able to provide detailed information on the molecular processes that underlie the technological application of materials in energy conversion and storage devices. Recent developments offer an increasingly predictive role for the techniques in material science generally, and energy materials in particular.
We thank our collaborators and colleagues for their insightful discussions, in particular, K. F. Aguey-Zinsou, A. A. Al-Sunaidi, J. S. Braithwaite, S. T. Bromley, S. Chen, S. F. J. Cox, S. Cristol, H. J. J. van Dam, D. C. Sayle, T.-X. Sayle, S. A. French, J. D. Gale, J. L. Gavartin, X. G. Gong, S. Hamad, M. S. Islam, E. Kaxiras, J. Kendrick, K. R. Kganyago, F. King, Y. Lei, N. M. Netsianda, P. E. Ngoepe, M. Nolan, S. C. Parker, M. Paul, S. C. Rogers, C. Schön, P. Sherwood, B. Slater, A. M. Stoneham, G. W. Watson, K. C. Waugh, S.-H. Wei, L. Whitmore, K. Wright, X. Xia, Z. Zhang and W. Zhu. We gratefully acknowledge funding by EPSRC (grants GR/T11364/01, GR/S26965/01, EP/E040071/1, GR/S52636/01, EP/E046193/1, EP/D504872 and EP/F067496) and Accelrys, and our membership of the Sustainable Hydrogen Energy and HPC Materials Chemistry Consortia. A.W. would like to acknowledge funding from a Marie-Curie Intra-European Fellowship from the European Union under the Seventh Framework Programme.
One contribution of 13 to a Discussion Meeting Issue ‘Energy materials to combat climate change’.
- © 2010 The Royal Society