Energy landscapes of perfect and defective solids: from structure prediction to ion conduction

The energy landscape concept is increasingly valuable in understanding and unifying the structural, thermodynamic and dynamic properties of inorganic solids. We present a range of examples which include (i) structure prediction of new bulk phases including carbon nitrides, phosphorus carbides, LiMgF3 and low-density, ultra-flexible polymorphs of B2O3, (ii) prediction of graphene and related forms of ZnO, ZnS and other compounds which crystallise in the bulk with the wurtzite structure, (iii) solid solutions, (iv) understanding grossly non-stoichiometric oxides including the superionic phases of δ-Bi2O3 and BIMEVOX and the consequences for the mechanisms of ion transport in these fast ion conductors. In general, examination of the energy landscapes of disordered materials highlights the importance of local structural environments, rather than sole consideration of the average structure.


Introduction
The underlying theme of this review is the energy landscape [1] and its use in computational solid-state science. This is an energy surface for the system of interest, a multidimensional plot or map showing how the energy depends on the variables that define the structure. Most usually these variables are the unit cell parameters and basis atom positions, sometimes referred to as external and internal degrees of freedom, respectively. The energy plotted is the internal energy (E or U) calculated using quantum theory within the Born-Oppenheimer approximation or molecular mechanics methods. Figure 1a, b shows schematic plots of a typical landscape. The traditional name in chemistry textbooks for such a landscape is "potential energy surface". The word "potential" is used as it is the form of the surface or landscape that defines the effective potential and this in turn determines the nuclear motions.
The lowest point in the deepest valley on the map corresponds to the global minimum, the lowest energy stable structure, while other valleys contain local minima corresponding to metastable structures with higher energies. Using this simple landscape analogy, it is straightforward to see that the dynamics of the system also relates to the energy landscape. The steepness of the valleys is related to the vibrational frequencies. The saddle points separating the valleys are transition states and indicate the lowest energy pathways between different possible structures. The energy barriers between the minima are also important; a structure corresponding to a local minimum will be metastable if the energy barrier between it and a lower energy structure is sufficiently high and thermally inaccessible. In disordered systems or solid solutions, many local minima may be thermally accessible, and understanding the properties of such systems requires thermodynamic averaging over these minima. At high applied hydrostatic pressures Published as part of the special collection of articles "Festschrift in honour of Prof. Ramon Carbó-Dorca". P ext , at which higher density structures are stabilised relative to those with lower densities, the enthalpy (H = U + P ext V) rather than the internal energy must be considered (Fig. 1c). Minima corresponding to stable structures at lower pressures may become saddle points in the energy landscape at higher pressures. Similarly, energy landscapes can be calculated for systems under uniaxial stress [2].
Further important considerations in structure prediction include temperature and, for systems which are not periodic in three dimensions, interfacial energies. At high temperatures, thermodynamic stability of different phases is determined by relative free energies as entropy cannot be neglected; this provides further substantial computational challenges, as the relative energies of different minima change. The energy landscapes of nanoclusters and ultrathin films can be very different from those of the corresponding bulk systems, since inclusion of interface energies can radically alter the shape and form of the energy landscape. Examples of these are given below.
Computationally, an enormous challenge is the large number of variables that need to be considered in order to achieve at best a very restricted exploration of the energy landscape, even for what might be thought of as "simple" systems with only a small number of atoms in the primitive unit cell. Ab initio finite-temperature free energy calculations are currently still very expensive. The problem of how to consider the topology of energy landscapes has been considered in detail by, e.g. Oganov and Valle [3]. In the energy landscapes of chemical systems, they note that low-energy minima are often found clustered together in a very few small regions or "funnels"-sometimes just one. Structural diversity is reviewed primarily in the context of phase transitions by Schön, Jansen and Doll [4].
In this short article, we cannot offer a comprehensive survey of the enormous amount either of the different methodologies or applications in this field, which are reviewed elsewhere (e.g. refs. [5][6][7]). Instead, we concentrate on a few representative examples taken from different areas and from our own work, presented in order of complexity. We start with the predicted energy landscapes of crystalline compounds, with most emphasis on binary compounds containing carbon. We then move on to the energy landscape of thin nanofilms of zinc oxides to show how the landscape of a system periodic only in two dimensions can be very different from that of the bulk material. We end by considering solid solutions and disordered systems, examining in particular two compounds containing Bi (δ-Bi 2 O 3 and BIMEVOX) which are excellent oxygen ion conductors. The examples discussed here illustrate the power of chemical intuition and "similarity" concepts in this context, as well as the importance of local bonding environments. We hope that our examples will be sufficient not only to illustrate the basic principles and insights obtained, but also to raise some more general issues.

Structure prediction: bulk materials
Often quoted and still highly apposite remarks, even after twenty-five years, are those of Maddox [8]: "One of the continuing scandals in the physical sciences is that it remains in general impossible to predict the structure of even the simplest crystalline solids from a knowledge of their chemical composition….". In the theoretician's defence, energy landscapes are highly multi-dimensional, and there is an exponential increase of the number of local minima with system size. Energy differences between different minima are often small and potentially unreliable given currently computationally feasible techniques for calculation of the energy [9]. A crystal structure is specified by the positions of all the atoms within a unit cell and by its lattice. Strain can alter both the shape of the lattice (specified by external coordinates) and the atomic positions (specified by internal coordinates). In order to find the most stable state of a crystal structure, i.e. to optimise a crystal structure, the appropriate thermodynamic potential is minimised with respect to a set of structure parameters. For more details and application to structure minimisation at finite temperature and applied pressure, see, for example, ref. [10].
For small unit cells, it is now feasible to use ab initio random structure searching to explore the energy landscape varying both unit cell parameters and basis atom positions [11,12]. This removes any danger of biasing the search towards structures based on known structures, but it is currently computationally unfeasible for larger systems. There can be a severe limitation on the number of atoms in the primitive unit cell. An illustrative example is an early application to N 2 at high pressures, for which a P − 4 2 1 m structure was predicted over a pressure range of approximately 10 to 56 GPa [12] with 8 atoms in the unit cell. Subsequent work of Ma et al. [13] illustrates the drawback of using such small primitive cells since a larger unit cell of 16 atoms and use of a genetic algorithm lead to a similar but not identical structure with Pba2 symmetry.
More generally, there are a number of diverse approaches to exploration of the energy landscape, as outlined, for example, by Woodley and Catlow [5], Oganov et al. [6], and Woodley et al. [14]. Methods include: • Simulated annealing [15] Here, a Monte Carlo (MC) algorithm together with the Metropolis criterion or short molecular dynamics runs generate a sequence of configurations across the energy landscape. Simulations start at high temperatures so that energy barriers between minima in different regions of the energy landscape can be crossed, and gradual cooling is carried out similarly to that in the physical process of annealing. As the temperature is lowered, low-energy regions are located. A related method uses the threshold algorithm [16]; • Basin hopping also involving passing through a sequence of configurations across the landscape but does so with the accept or reject decision made on the magnitude of the energies of the local minima visited. Each trial move in the landscape will involve a full structural relaxation following either a Monte Carlo step or a short molecular dynamics run [17]; • Evolutionary (genetic) algorithms in which a population of structures evolves as lower-energy structures are generated by natural selection from an initial set of structures by mutations, genetic crossover and replication [18]; • Metadynamics in which all structural variables are reduced to a collective variable capable of distinguishing different structures. Then, low-energy structures and pathways and the energy barriers between them are found; as it proceeds, the simulation is discouraged from sampling regions of the landscape already sampled [19]; • Topological methods for exploring networks and particularly useful for framework structures such as zeolites [20] and metal-organic frameworks (MOFs) [21].
Methods based on machine learning are rapidly becoming more important [14]. All methods rely on the accurate and rapid calculation of energies either by quantum-mechanical methods, chiefly DFT, or by interatomic potentials (force fields), which are increasingly developed by machine learning trained on quantum-mechanical results for smaller systems.
An alternative to global searching methods is identification of possible candidate structures by data mining followed by refinement by quantum mechanical structure optimisation of both unit cell parameters and basis atom positions. As well as data mining, intuition can often be useful in identifying candidate structures and chemical trends which can indicate further possibilities [22]. In a related approach, known structures for one of the constituent elements in its pure form can be taken as a starting point and modified to include the other elements. While not strictly a first principles method, this is a promising way of probing interesting parts of the energy landscape given the conclusions of Oganov and Valle [3] regarding the typical form of energy landscapes and the "funnels" therein. For example, in our work on carbon-containing compounds, we have often used diamond and graphite as starting points and systematically introduced nitrogen or phosphorus, making sure that chemical valency rules are satisfied. Kroll and Hoffman [23] have considered possible structures for Si 3 B 3 N 7 in this way, deriving initial structures by systematic introduction of boron into the structures known for α-and β-Si 3 N 4 . It is important to stress that in our work, we use large supercells containing multiple formula units, and both unit cell parameters and the positions of all basis atoms are allowed to vary during structural geometry optimisations so that even if the initial guess involves a known structure, the final minimised structure can be very different from the starting candidate structure and a completely new crystal structure thus generated. The data mining and chemical intuition serve simply to find realistic starting structures with appropriate stoichiometries, distances, and chemically sensible local coordinations and environments.
To turn to some examples, there has long been considerable interest in carbon nitrides [24][25][26][27]. Some binary carbon-nitrogen compounds have been predicted to be harder than diamond, and others have potential applications as photocatalysts [28] and as semiconductors [29]. The lowest energy structures at ambient pressure are all graphiterelated. For CN and C 3 N 4 , the preferred structures contain six-membered rings and are related to graphite, with some carbon replaced with nitrogen and vacancies ensuring that all valencies are satisfied. Structures can be either twodimensional sheets or one-dimensional chains; some low, virtually equienergetic examples for CN found in our work are shown in Fig. 2a, b. For C 3 N, the lowest energy structure found by Hart et al. [22] and by Sandré et al. [30] is similarly obtained by starting from graphite. In contrast to CN, however, valencies are satisfied by the formation of C-C bonds between layers (Fig. 2c), rather than the introduction of vacancies. After optimisation, we found that there are considerable deviations from planarity for some of the atoms; some six-membered rings are close to planar, while others are more chair-like. For CN at higher pressure (and incidentally also for SiN and GeN at ambient pressure), we found a structure derived from the β-InS structure to be lowest in energy (Fig. 2d) [22]. The β-InS-like structure consists of tetrahedral C/Si/Ge bonded to one C/Si/Ge atom and three almost planar threefold coordinate nitrogens. For any given carbon nitride composition, there are several possible structures very close in energy, often with commonalities in local bonding environments; nevertheless, there can be appreciable differences in density (e.g. compare Fig. 2a, b). For a perspective of the current experimental position, see Miller et al. [31].
Less attention has been paid to phosphorus carbide, also commonly referred to as carbon phosphide even though phosphorus is the more electropositive of the two elements. The low-energy structures at different P:C ratios are interesting contrasts to those in the C-N system, reflecting the different chemistries and structural preferences of phosphorus and nitrogen. There are significant changes with composition in the preferred local environments for both carbon and phosphorus [32,33]. For carbon-rich compositions (50 at% carbon), structures derived from graphite with four-coordinate, hypervalent phosphorus and sp 2 carbon are preferred (Fig. 2e). For higher phosphorus content (56 at% phosphorus), defect zinc-blende structures derived from diamond, containing three-coordinate phosphorus and sp 3 carbon, are lowest in energy (Fig. 2f). An experimental study of the synthesis of phosphorus carbide by pulsed laser deposition reports disordered films and broadly confirms the predicted local bonding environments [34]. As for carbon nitride, several possible structures are often very close in energy for a given composition, which is reflected in the difficulty of preparing crystalline forms of these and related systems. But we stress that for a given composition, the local environments of the carbon and phosphorus in the lowest energy structures are the same. Recent interest in this system has moved to carbon phosphide monolayers for nanoelectronics applications such as gas sensors [35].
We stress again that it is important to explore the energy landscape around the minima corresponding to known structures. This may well locate further minima which themselves do not correspond to any structure in the initial set or indeed to any known structure. Full geometry optimisations without symmetry constraints are crucial in increasing the chance of finding these new minima as significant bond making and breaking and structural changes can occur during optimisation. Optimisations initially under applied positive or negative hydrostatic pressures can also be useful, since this may allow kinetic barriers to be overcome and transitions to new structures not in the initial set of starting structures may take place. For example, the large deviation from planarity for six-membered rings in the lowest energy C 3 N structure was only seen after optimisation first under a negative hydrostatic pressure and then reoptimisation at ambient pressure [25]. In addition, once preferred local environments are established, it is often straightforward to see other possibilities. Figure 2g shows our predicted lowest energy structure for a ternary compound-LiMgF 3 [36]. Calculation also indicates a negative enthalpy of formation from the binary fluorides LiF and MgF 2 . The structure predicted is that adopted by LiNbO 3 and it is ferroelectric. Whereas this particular Fig. 2 Predicted low-energy structures for a 2D and b 1D forms of CN; c C 3 N (carbon grey, nitrogen blue); d SiN, GeN and high-pressure CN (Si/Ge/C grey, nitrogen blue); e PC; f P 4 C 3 (C grey, P orange); g LiMgF 3 (Li green, Mg purple, F red) structure is common for ternary oxides and even predicted by simple radius ratio arguments, it is intriguing that no fluoride with this structure has yet been reported experimentally.
There are a number of strategies such as high pressure for exploring experimentally different parts of energy landscapes. As mentioned earlier, at such pressures, the enthalpy must be evaluated, and structures with smaller volumes and thus higher densities are stabilised. The classic textbook example is the adoption by NaCl and other alkali halides of the CsCl structure in which the ions are eightfold coordinate at high pressures, above approximately 27 GPa for NaCl [37]. Textbooks make less mention of stable structures with unexpected stoichiometries formed by sodium and chlorine at higher pressures such as Na 3 Cl, Na 2 Cl, Na 3 Cl 2 , NaCl 3 and NaCl 2 , predicted and confirmed by diamond anvil cell experiments [38]. In recent years, dramatic advances have been made in predicting, and then preparing experimentally, high-pressure structures including, for example, new high-temperature superconductors [39] and the unexpected electride Na 2 He [40]. Another highlight has been the discovery of a partially ionic high-pressure phase of elemental boron through a combination of high-pressure experiments and evolutionary structure prediction [41]. Other examples include successful synthesis of Na 3 N [42] after a series of computational studies using simulated annealing and the threshold algorithm [43] over a wide pressure range. The inability to synthesise any compound with this composition had received considerable attention for many years because of the ease of preparation of the homologue Li 3 N. It is gratifying that almost the complete set of the most stable polymorphs of Na 3 N, predicted as a function of pressure, have now been prepared successfully [44].
Somewhat less attention has been paid to predictions of low-density polymorphs. Nevertheless, attention has been drawn computationally to possible nanoporous structures with lower densities than those of phases well-known under ambient conditions [45][46][47][48]. For example, ab initio calculations support the existence of new low-density nanoporous crystalline phases of MgO and ZnO which may be experimentally accessible by coalescence of nanocluster building blocks. Postulated structures have unit cell volumes 20-110% larger than those of known structures. Nanoporous phases of alkali metal halides have also been proposed [49].
Possible nanoporous low-density forms of B 2 O 3 containing B 3 O 3 rings have been linked to the ease of vitrification of this system [50] and found to be extremely flexible [51,52]. Starting with some known B 2 O 3 structures, and others constructed from two-dimensional sheets, we explored the effect of replacing each three-coordinated B atom in these polymorphs by a B 3 O 3 ring, itself three-coordinate. This generates new superstructures with the same B 2 O 3 stoichiometry, in which the rigid boron-oxygen B 3 O 3 heterocycles are linked by flexible B-O-B bridges. The basins in the energy landscapes of these framework structures are extremely shallow, with widths as large as those of many hybrid metal-organic frameworks (MOFs). Very large changes in cell volume-more than 200%-result in only negligible changes in energy, and the energies are comparable to those of known phases of B 2 O 3 . This flexibility is unprecedented for systems which are purely inorganic; we suspect many of these structures are present in the vitreous oxide given the extensive controversy over its density [53]. This is in marked contrast to what is known or has been predicted for low-density frameworks of the other materials mentioned above, where the lower the density the more destabilised is the structure relative to a denser ground-state polymorph.
In an extension of this work, similarity arguments indicate that given any arrangement of SiO 4 tetrahedra in a solid silicate, it is possible to construct a new structure in which the SiO 4 4− tetrahedra are replaced by B 2 O 5 4− units and in which all the boron atoms are three-coordinate. Vast arrays of theoretical B 2 O 3 polymorphs may thus be generated from known silicate compounds, including new families of zeolites, and we anticipate a general increase in flexibility.

Ultra-thin films: graphene forms of ZnO and other wurtzites
The intense attention paid in recent years to thin films and nanoclusters is revealing a wealth of new structures and associated chemical and physical behaviour, and here we outline briefly an example which illustrates strikingly the differences between the energy landscapes for bulk and thin films. Many II-VI and III-V semiconductors adopt bulk structures in which all the atoms are tetrahedrally coordinated, such as wurtzite or zinc-blende. These are wide band-gap materials with a wealth of potential applications in optoelectronic devices. The energy landscape of bulk ZnO in particular has been studied extensively. Under ambient conditions, the wurtzite structure (Fig. 3a) is adopted; the zinc-blende structure is metastable (Fig. 3b). At 9 GPa [54], ZnO undergoes a phase transition to the NaCl structure. As mentioned above, ab initio calculations at negative pressures have also discovered energy minima corresponding to more open structures based on zeolite motifs [45]. Ultra-thin films and nanoscale synthesis also offer important routes by which different structures can be stabilised. We have reported density functional theory (DFT) [55,56] and more recently hybrid B3LYP calculations [57] that demonstrate that a single layer of ZnO adopts a graphene-like structure, while for a small number of layers, an "eclipsed" graphite-like structure, isostructural with hexagonal boron nitride, is adopted in which Zn atoms in one layer lie above the O in that beneath and vice versa (Fig. 3). Substitution of all the atoms in the zinc-blende and wurtzite structures by carbon generates the face-centred cubic and hexagonal Lonsdaleite diamond structures, respectively; by analogy, the formation of a graphene form of ZnO might then be anticipated. The graphitic form in the bulk is unstable with respect to the wurtzite structure by approximately 0.3 eV per formula unit. Finite size effects, favouring the removal of the dipole moment perpendicular to the layers in the wurtzitestructured film, alter the relative stabilities for thin films. For discussion of the effects more generally of the polarity of oxide surfaces and nanostructures on energy landscapes, see the work of Noguera and co-workers (e.g. ref. [58]).
We have suggested that growth of ZnO crystals, the morphology of which is dominated by the {0001} surface, can begin with thin films with this "eclipsed" graphitic structure that then undergo a barrierless transition to the wurtzite structure when a critical number of layers is reached [55]. Subsequent to the original prediction, the existence of such graphitic structures was confirmed experimentally using electron microscopy [59].
Subsequent studies have shown that the structural chemistry of ultra-thin films of ZnO is even richer. Other lowenergy possibilities for a small number of layers, shown also in Fig. 3, include a body-centred tetragonal (BCT) form [60] comprised of Zn 3 O 3 rings which are not flat but have boat conformations and a "staggered" graphitic form [57] in which adjacent graphene sheets are displaced relative to each other. The latter is particularly interesting in showing a spontaneous symmetry breaking across the film and a band gap that decreases substantially with film thickness, associated with the transfer of electron density from one layer to another.
Other compounds with the wurtzite structure should also demonstrate this behaviour [56]. For ZnS, which is most stable in the zinc-blende structure, the low-energy polymorphism is rather less extensive than for ZnO, related to the relative preference of S compared to O for fourfold coordination over threefold planar coordination; polar ZnS nanofilms are predicted to lose Zn ions spontaneously and so become non-stoichiometric [57].

Energy landscapes of disordered materials
For a periodic solid such as ZnO, different minima in the bulk potential energy landscape correspond to different possible crystal structures. In addition, for a solid solution, a disordered solid or a non-stoichiometric compound, the energy landscape concept can also be highly instructive and the remainder of this article is devoted to some typical examples.

Solid solutions
In a substitutional solid solution such as that between binary halides [61] or binary oxides [62,63], the energy of each possible arrangement of the cations is associated with a Fig. 3 Possible structures of ZnO films: (derived from) bulk wurtzite; (derived from) bulk zinc-blende; graphitic (graphene sheets "eclipsed" such that ions in one sheet lie directly above ions of opposite charge in the sheet below); alternative graphitic form "staggered" (graphene sheets "staggered" such that adjacent sheets are displaced relative to each other); BCT (body-centred tetragonal). Zn ions blue, O. "Eclipsed" graphitic is the lowest in energy for ultra-thin films different local minimum, and we can build up an energy landscape showing these energy minima and the corresponding atomic arrangements. It is often merely assumed that the different atoms are distributed at random over the available sites (the ideal solution model); this corresponds to the energy landscape in Fig. 4a where all the minima have the same energy independent of the detailed arrangement of the ions. Often this is not so, such as when the different atoms vary considerably in size, and solutions can be strongly nonideal, and the atoms are far from randomly distributed; in this case, different arrangements have different energies as in Fig. 4b and Fig. 4c. By analogy to geographical landscapes, one can imagine many possibilities-some rugged and complex, others smooth and simpler in form. The type of disorder is strongly connected to the type of energy landscape. To calculate the thermodynamic properties of the solutions to assess the empirical assumptions often made (e.g. regular solution theory, assumption of ideal entropies of mixing), we need to carry out thermodynamic averaging over the different arrangements (configurational averaging). In principle, the solid solution can assume any state in which each atom can be at any position. However, the only states of practical importance away from the melting point will lie at the local minima at the bottoms of the basins in the energy landscape, each corresponding to a given configuration and we need to sample only these basins. Using the label k = 1…K for the configuration, then the enthalpy, H, and Gibbs energy, G, in the isobaric-isothermal (NPT) ensemble are given [63] by where k B is Boltzmann's constant. H k and G k are the enthalpy and Gibbs energy, respectively, for the relaxed structure of each possible cation arrangement. We thus have expressions for any thermodynamic quantity in terms of thermodynamic quantities obtained for particular configurations. In favourable cases [62,63], the thermodynamic averaging can even be performed over the results of a set of full free-energy minimisations [64] of different arrangements (configurations) of the cations within a supercell. To calculate G k , the quasiharmonic approximation can be used to calculate the phonon frequencies over sufficient wavevectors in the Brillouin zone to ensure convergence. In practice, G k is often replaced by H k in Eq. (1) and Eq. (2), and the vibrational entropy and all vibrational contributions to the enthalpy are neglected.
A serious difficulty with the practical use of these equations is the large value of K., the number of minima in the energy landscape. Other than with the smallest cells, which may exclude the most crucial arrangements, K is impractically large. For intermediate cell sizes, use of symmetry via explicit identification of the space group of each supercell arrangement [63] is highly valuable for finding the weighting of each state of each minimum, but in many cases, we must restrict the sum to an incomplete but representative sample of minima. In strongly non-ideal situations, where only a few deep local minima in the energy landscape are thermally accessible, a genetic algorithm is often useful [65][66][67] to identify the low-lying orderings [68]. We have also developed a set of Monte Carlo techniques [63,69] for these problems, which involve explicit exchanges of atoms and steps designed to allow for structural relaxation of the individual configurations.
We refer the reader to the references for full details of the configurational averaging, but it is worth highlighting some of the conclusions. Formally, this is a "basin sampling" technique [70] which has been used for biomolecules [71,72] and clusters [73] as well as for solids. The local minima are often called "inherent structures" when considering glasses and supercooled liquids [74][75][76].
For solid solutions of halides or oxides, the solutions are non-ideal unless there is only a very small size mismatch between the cations on the same sublattice; clustering will occur to some extent. Consequently, the configurational entropy of mixing can be substantially less than the ideal value ( R ∑ i x i lnx i ), which is the value typically assumed. But the total entropy of mixing-configurational plus vibrational-can be considerably higher than the "ideal" value. For CaO-MgO at 2000 K, the total entropy of mixing at 50-50 composition is 50% higher than the "ideal" [69]. Examination of the structures corresponding to local minima in the energy landscapes is also instructive allowing local structure to be related to energetics. For example, in pyrope-grossular (Mg 3 Al 2 (SiO 4 ) 3 -Ca 3 Al 2 (SiO 4 ) 3 ) garnet solid solutions, low-lying minima in the landscape correspond to specific orderings of the Ca 2+ and Mg 2+ ions. Surprisingly third and fourth nearest neighbour cation interactions dominate the orderings and in each case Mg-Mg pair combinations are energetically unfavourable. For the third nearest neighbours, this has been linked to local distortions of the SiO 4 tetrahedra by the smaller cation [77] and for the fourth neighbour interactions to distortions of AlO 6 octahedra. In other cases, such subtle changes in local geometry are important more generally in determining the sign of the enthalpy of mixing of two solids and hence whether a solid solution or an ordered phase will form. An interesting contrast is that between the oxides CaO-MgO, where the enthalpy of mixing is highly endothermic for all compositions, and the corresponding carbonates CaCO 3 -MgCO 3 which despite the similarity in crystal structure combine exothermically to form the ordered dolomite CaMg(CO 3 ) 2 . Examination of the structures corresponding to different minima in the energy landscape [78], located by a genetic algorithm, shows that this very different behaviour is thanks to relief of strain in the Mg-O and Ca-O bonds made possible by suitable rotations of the polyatomic carbonate anions.

Highly non-stoichiometric compounds
At first sight, non-stoichiometric compounds containing vacancies are a close relative of the substitutional solid solutions. Instead of distributing different atoms on a given sublattice, one type of atom and vacancies are distributed on a sub-lattice. While this can serve as a useful first model, in practice for oxides, we have found these systems are structurally much more closely related to many amorphous systems in that a strong preference for certain coordination polyhedra dictates the local order in the overall disordered state. Structural relaxations accompanying different distributions of these polyhedra may be so large (as in Ba 2 In 2 O 5 ) [79] that the idea that vacancies occupy a sub-lattice is a poor model of the disorder although this is commonly accepted practice. It can instead be useful to compare the picture of non-stoichiometric compounds with that of an amorphous material such as vitreous silica. In such cases, the atoms are not distributed on a well-defined sub-lattice, but there is considerable short-range order, in the form of SiO 4 tetrahedra, and intermediate-range order arising from the interactions between the polyhedra, revealed by such indices as the Steinhardt order parameter [80]. This type of model is related to the theory of supercooled liquids and the landmark paper of Goldstein [81]. In undercooled liquids, the atoms will spend most of their time in an arrangement corresponding to a local minimum in the energy landscape. Occasionally, a transition to a neighbouring basin will take place. At higher temperatures intervalley jumps are more common, but are still rare compared to intravalley movements. The methods of the previous section and Eqs. (1) and (2) can be used to extract the thermodynamic properties. Similar such methods have been used, for example, to probe the relationship between fragility, configurational entropy and the energy landscape of glass-forming liquids [82]. Eventually at yet higher temperatures, the separation of the partition function into separate configurational and vibrational parts will no longer be valid.
For grossly non-stoichiometric compounds, solid solutions or disordered materials more generally, most crystallographic techniques produce an 'average' structure based on partial occupancy of particular sites (see the schematic plot in Fig. 5). These yield no information about which local atomic arrangements or configurations give rise to this average. In the absence of more information and in simple statistical mechanical treatments such as mean-field theory or regular solution models, it is usual to assume random distributions of cations or vacancies consistent with this average. But such assumptions can fail dramatically.
For example, studies of the energy landscapes of oxygendeficient perovskites such as Sr 2 Fe 2 O 5 and Ba 2 In 2 O 5 reveal a strong energetic preference for just a few local structural entities, including Fe(In)O 4 tetrahedra, Fe(In)O 6 octahedra, and some Fe(In)O 5 pyramids [79,83]. The average structure observed experimentally at elevated temperature can then be interpreted as a spatial average of these local structures involving tetrahedra, octahedra and pyramids.

Fast-ion conduction -LAMOX, Bi 2 O 3 and BIMEVOX
The importance of local structure has major implications for understanding both thermodynamic properties and dynamics. Ion transport in the fast-ion oxygen-deficient perovskites, for example, involves transitions between the basins corresponding to different local minima in the energy landscape and thus the interplay of these structural entities. Thus, we explore via their energy landscapes the links between local structure and observed physical and materials behaviour, including structure, thermodynamics, and transport. This project is challenging due to the difficulties associated with the modelling of solid solutions and grossly non-stoichiometric systems once one moves away from assuming ideal or regular solutions or defect distributions or from traditional textbook 'mean-field' treatments which fail to allow for distinct local environments of the different species present. The challenge computationally is that of calculating the energies of a large number of possible atomic arrangements while taking full account of the changes in local environment of each species from one arrangement to another. There arises an inevitable balance between (computational) expense and the number of arrangements considered.
Using the same general approach, correlated movements of the oxygen atoms in the important superionic conductor LAMOX, La 2 Mo 2 O 9 , have been demonstrated by Hou et al. [84] and these highlight the significance in LAMOX of certain four-or fivefold coordinated molybdenum. The importance of local atomic arrangements in LAMOX and the implications for the monoclinic-to-cubic phase transition have also been elegantly probed experimentally by Malavasi Fig. 5 Schematic picture of two local structures. The "average" structure is in fact not a local minimum but a transition state et al. [85] Cooperative mechanisms of fast-ion conduction have also been established [86] in gallium-based oxides with tetrahedral moieties such as LaBaGaO 4 , where a cooperative "cog-wheel" process operates, in which Ga 2 O 7 units form and break up. Correlated mechanisms are emerging as much more important in fast-ion conduction than previously recognised, even in such well-studied superionics such as PbF 2 where previously there was little evidence for collective diffusion [87].
Oxide ion conductivity has been investigated much more thoroughly than hydride ion conductivity, a much more recently observed phenomenon that has so far been restricted to a small number of recent studies of metal hydrides and oxyhydrides. We have recently studied hydride transport in Ba 2 ScHO 3 , a recently synthesised oxyhydride with an unusual anion ordering, examining with DFT the minima in the energy landscape and the transition states between them and also carrying out ab initio molecular dynamics simulations [88]. Ionic mobility depends on the hydride-oxide disorder but activation energies to migration do not; understanding how the migration energies depend on local structural flexibility has allowed us to predict that Sr 2 ScHO 3 will have larger conductivity than the Ba analogue. Through the study of energy landscapes, we are thus now able to develop principles of crystal engineering not only for new structures, but also for specific dynamic properties such as ion conduction.
We turn now to consider some fast ion Bi compounds, which are highly non-stoichiometric and of considerable current interest. We start with the disordered δ-phase of bismuth oxide, Bi 2 O 3 , which has very high ionic conductivity. While the superionic phase is only stable at high temperatures (approximately 1000-1100 K), it can be stabilised even at ambient temperatures on SrTiO 3 or DyScO 3 [89]. The thermodynamically stable bulk phases of bismuth oxide (α-and β-) at lower temperatures are ordered and have much lower conductivity. Sr 2 Fe 2 O 5 and Ba 2 In 2 O 5 can be viewed as oxygen-deficient perovskites; δ-Bi 2 O 3 similarly can be described formally as a Bi 4 O 8 fluorite lattice with one-quarter oxide sites vacant. However, as for the pervoskites, this simplification can be highly deceptive. Experimentally, diffuse neutron-scattering combined with DFT configurational averaging and Born-Oppenheimer ab initio molecular dynamics [90][91][92] reveal a very different scenario from the conventional picture of an average cubic structure in which oxygen vacancies are distributed randomly. Rather, the local Bi environments in δ-Bi 2 O 3 are very similar to the distorted square pyramids present in the ordered β-phase. In the δ-phase, there are many different orientations of these pyramids and there are many local minima in the landscape corresponding to these different orientations of distorted Bi environments, separated by low-lying barriers, as shown schematically in Fig. 6a. Only the average structure is cubic; transitions between the energetically accessible minima in the rugged energy landscape, which correspond to low-symmetry structures, are similar to those in the ordered phases. These are responsible for facile oxygen conduction in δ-Bi 2 O 3 . In contrast, the "average" cubic structures, also shown in Fig. 6a, are all high in energy and thermally inaccessible. In disagreement with experiment, these cubic structures are all metallic or semi-metallic at the GGA level of theory, while the lowenergy structures are all insulators with a calculated GGA band gap of ~ 2.0 eV which compares with an experimental value of 2.6 eV.
Consideration of the bond lengths alone is sufficient to illustrate the striking difference between the local coordination environments in δ-Bi 2 O 3 and the mean structure. Average Bi-O and Bi-O / distances in this phase are ~ 2.45 Å and ~ 2.33 Å, respectively [91]. But in the ordered, low temperature α-and β-phases [91] containing the same Bi local environments and coordinations, only 10% and 20%, respectively, of the Bi-O distances lie in the range 2.40-2.50 Å. There are no Bi-O separations at all in the range 2.30-2.40 Å. Overall, there are very few bonds indeed in δ-Bi 2 O 3 with the mean bond length. These variations in actual bond lengths from the average and the distorted local environments reflect the effect of the asymmetric lone pair on Bi. It is tempting to associate the facile ion conduction with this variation and the ease of distortion of the local environment. Some support for this comes from examining the energy landscape and corresponding local Bi environments of Bi 2 O 3 adsorbed at the (100) surface of SrTiO 3 compared with those for the bulk [93]. Low temperature stabilisation of the superionic phase is not due to any epitaxial matching, but rather the opposite, i.e. the mismatch between Bi-O and Sr-O bond lengths. Small Bi-O domains at the SrTiO 3 /Bi 2 O 3 interface promote the formation of a disordered layer of δ-Bi 2 O 3 , with irregular orientations of the same local environments of the Bi as those in the ordered phases. There is no stabilisation of a symmetric cubic structure. BIMEVOX, i.e. γ-Bi 4 V 2 O 11, and its derivatives have exceptionally high ionic conductivity even at 500 K [94]. Two lower temperature phases (α-and β-) have much lower conductivities. In the γ-phase, Bi 2 O 2 2+ layers are separated by perovskite-related VO 3.5 2− slabs; the structure can be generated by the introduction of oxygen vacancies into the "parent" Bi 2 VO 6 Aurivillius structure; we used a 2 × 2 × 1 supercell of the "parent" introducing eight O vacancies. Different oxygen vacancy arrangements are shown in Fig. 6b together with the key features of the corresponding energy landscape calculated using DFT [95]. In this figure, configurations are grouped according to the location of the oxygen vacancies; the labels 't' and 'b' refer to vacancies being in adjacent vanadium layers and 'm' and 'e' to vacancies in Once more, there are many minima in the energy landscape of γ-Bi 4 V 2 O 11 , dominated by low-lying thermally accessible 4t4b configurations; examination of the optimised geometries shows a corresponding range of vanadium coordinations and large variation in the Bi-O and V--O bond lengths. The grouping of local minima in the energy landscape into sets of configurations also permits us to examine diffusion separately in the different layers using ab initio molecular dynamics simulations and then to calculate properties such as the average diffusion coefficient by configurational averaging analogous to that for the enthalpy in Eq. (1). Oxygen transport is cooperative, involving transport of oxygen ions in curved paths through the vanadium layers [95], confirming experimental evidence that oxygen disorder is confined to these layers [96]. The conduction mechanism is thus very different in BIMEVOX from that in δ-Bi 2 O 3 . In δ-Bi 2 O 3 , Bi ions can be coordinated to 4, 5, 6 or 7 oxygens and transport involves transfer of oxide ions bonded to Bi and changes in the bismuth coordination. On the other hand, in Bi 4 V 2 O 11, it is the oxygens bonded to vanadium that have high mobility and the coordination numbers of the vanadium vary more than those of the bismuth. That said, oxide transport is helped by the ease of distortion of the local environments of the Bi, shown by the range of Bi-O separations thermally accessible and the curved migration pathways. Once more, the importance for ion conduction of the local environments present in the disordered phase is demonstrated.

Electronic Properties
The ability to study the thermodynamics and structures of solid solutions and disordered structures prompts the question as to what insights can be gained into electronic properties and how these may vary or be optimised for particular applications. For example, it is often highly desirable to be able to tune the band gap of semiconductors. Many photocatalytic systems for the generation of hydrogen from water unfortunately have band gaps too large for the absorption of visible light and hence are not efficient for use under sunlight, where the most intense wavelengths are in the visible region. Thus, there is much interest in modifying band gaps to achieve photocatalysis in this region. One approach is via combining semiconductors to form solid solutions. There are many examples in the literature; in our own work, we have studied [97] the GaP-ZnS solid solution and how the band gap varies with atomic ordering and composition, with hybrid functionals (B3PW and HSE06). It is interesting to note once again the importance of local bonding environments, now in terms of electronic rather than structural effects -the energies of states at the top of the valence band (anion 3p-states) can be particularly sensitive to small changes in electrostatic potential and thus local environment, corresponding to significant changes in band gap. GaP itself has an indirect band gap (and so is inefficient for light absorption), while the band gap of ZnS, although direct, lies in the ultraviolet region. Configurational averaging of the band gap [98] (analogously to Eq. 1) suggests that for (ZnS) 1-x (GaP) x compositions where x ≤ 0.75, the band gap is direct and we predict that compositions with x = 0.25, 0.5 and 0.75 are promising photocatalysts for the production of hydrogen from water under solar radiation due to their band gaps and the positions of the bands relative to vacuum. The potential for achieving visible-light activity with the ZnS-GaP system due to control of local bonding environments has been confirmed in experimental work [99].

Conclusions
The energy landscape is an increasingly valuable concept for understanding the structural, thermodynamic and dynamic properties of inorganic solids, as it also has proved to be for small molecules and large biomolecules such as proteins and enzymes [100]. In this article, we have covered a range of examples, including structure prediction of new bulk phases and ultra-thin films to understanding highly nonstoichiometric oxides and solid solutions at high temperature. Increasing attention is being paid to the transition states (saddle points) in the landscape, which govern the kinetics, diffusion and conductivity and crystal growth.
In general, examination of the energy landscapes for both the bulk phases and the disordered systems highlights the local structural environments which are present. Thus, consideration of the average structure alone is not sufficient, and as we have seen can at worst be highly misleading. In solid solutions and non-stoichiometric compounds, non-ideality is strongly linked to the local environment, and our results question the use of models which assume completely random distributions of ions or vacancies over particular lattice Fig. 6 Schematic energy landscapes of a δ-Bi 2 O 3 (note especially the very high energy of the three "average" cubic structures; < 100 > , < 110 > and < 111 > denote different orientations of the two oxygen vacancies in the "parent" cubic fluorite Bi 2 O 4 ), and b BIMEVOX -α-, β-and γ-Bi 4 V 2 O 11 . In the inset showing the structure of the α-phase, Bi atoms are purple, V large red, O small red. The labels '4t4b', '6t2b', '3t3b1m1e' and '2t2b4m' refer to different oxygen vacancy arrangements: 't' and 'b' refer to adjacent V layers, 'm' and 'e' to adjacent Bi layers and the numbers indicate the relative numbers of vacancies in each layer. For example, 4t4b indicates there are equal numbers of vacancies in all the V layers and none in the Bi layers. The vanadium-oxygen polyhedra are shaded, and for clarity Bi, atoms are not shown ◂ sites. There are important consequences for transport and correlated ion motion as well as thermodynamic behaviour, electronic properties and solid-state reactions such as trace element (dopant) incorporation. With advances in computer power and ab initio computational methodologies, the long-term prospects of our work are a step change in our understanding of how the atomic-level behaviour governs the macroscopic physical and chemical behaviour, and of the importance of the local structure and environment in determining materials-related properties of technological interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.