Adaptive kinetic Monte Carlo simulation of solid oxide fuel cell components

Abstract

Ionic conductivities in the solid oxide fuel cell (SOFC) electrolytes yttria-stabilised zirconia (YSZ), calcia-stabilised zirconia (CSZ), gadolinium-doped ceria (GDC) and samarium-doped ceria (SDC) and the cathode material lanthanum strontium cobalt oxide (LSCO) are directly calculated using DL_AKMC, an adaptive kinetic Monte Carlo (aKMC) program which assumes limited a priori knowledge of the kinetics of systems. The materials were simulated over several milliseconds and over the range of experimentally most relevant temperatures and dopant concentrations (2–18 mol% for doped zirconia, 5–25 mol% for doped ceria and 5–80 mol% for LSCO). Ionic conductivities of the electrolytes at 1000 K are in good agreement with the observed values: CSZ in the range 3 × 10⁻³ to 1 × 10⁻² S cm⁻¹ depending on dopant concentration, YSZ 4 × 10⁻³ to 3 × 10⁻² S cm⁻¹, GDC 1 × 10⁻² to 5 × 10⁻² S cm⁻¹, SDC 1 × 10⁻² to 7 × 10⁻² S cm⁻¹. LSCO is predicted to have an ionic conductivity of the order of 10⁻² to 10⁻¹ S cm⁻¹ depending on Sr content. Average activation energies over all migration processes are 0.4–0.5 eV for the stabilised zirconias and 0.2–0.3 eV for the doped cerias and 0.3 eV for LSCO, in agreement with experiment. aKMC provides a distinct advantage over traditional KMC methods, in which one has to provide a list of system state transitions. Here, all of the state transitions are dynamically generated, leading to a more accurate simulation of the kinetics as the system evolves.

---

1 Introduction

Due to growing concern about global warming and dwindling supplies of fossil fuels, increasing interest is being directed towards alternative sources of energy. Popular and promising alternatives are fuel cells, and a substantial amount of research and development has been undertaken to improve the materials used in these devices. In solid oxide fuel cells (SOFCs) the oxidant (e.g. air, O₂) is reduced at the cathode. The oxide ions produced are then transported through the solid electrolyte material, ideally a purely ionic conductor, to the anode where the fuel (e.g. H₂, hydrocarbons) is oxidised. Electrons from this process then flow from the anode to the cathode, completing the circuit and generating power. Lowering the operating temperature of SOFCs to an intermediate temperature range of 600–800 °C is an ongoing area of intense research and many materials have been suggested as suitable electrolytes^1–3 and electrodes^4–9 for this purpose.

Stabilised zirconias doped with lower valency ions such as yttrium (YSZ) or calcium (CSZ) along with ceria doped with gadolinium (GDC) and samarium (SDC) exhibit very high ionic conductivities associated with mobile oxygen vacancies formed as a consequence of the doping and are typically used as electrolytes in SOFCs. Previous computational studies of YSZ have tended to use molecular dynamics (MD)¹⁰ or density functional theory (DFT)^11,12 in combination with kinetic Monte Carlo (KMC) to investigate the oxide ion diffusion kinetics,^13,14 as timescales tend to be limited if MD or particularly DFT is used in isolation. Typically common atomistic KMC models employ an on-lattice approximation which limits their ability to describe a system which undergoes large structural changes. This approach also requires a list of possible event mechanisms determined a priori, through experimental and theoretical methods, by estimation or even guessing. This is naturally a severe limitation as the processes involved in the atomic motion are not necessarily intuitive and can be extremely difficult to predict in advance.¹⁵ Furthermore, as the simulation advances the structure of the system changes leading to new possible transitions and altering the activation energies of existing transitions. To overcome these limitations various ‘on-the-fly’ approaches, such as adaptive KMC (aKMC),¹⁶ have been proposed. These are designed for an off-lattice system and the saddle-points between different possible states of the system are located as the simulation progresses. These methods overcome the limitations of the traditional KMC approaches by avoiding the need for a list of mechanisms provided beforehand. State transitions not predicted in advance are allowed, thus permitting a thorough search of the energy landscape and ensuring that the kinetic model is as realistic as possible. As in traditional KMC, certain assumptions are made when using aKMC, namely that the system transitions adhere to harmonic transition state theory.^17,18 Here, we also assume that the pre-exponential factors, used in the calculation of transition rates, are equal to a typical vibrational frequency (10¹³ s⁻¹).

2 Methodology

2.1 aKMC

The key to the aKMC approach is an efficient way to search for saddle points that link the current state of the system to another. One popular approach is the ‘dimer’ method,¹⁹ the implementation of which used in this work has been described in detail elsewhere.²⁰ The dimer algorithm begins with a random starting position within the energy basin and progresses by climbing uphill along the lowest eigenvector corresponding to the lowest value eigenvalue of the Hessian matrix, eventually reaching a saddle point at the top. This is a particularly efficient method as only first derivatives of the energy²¹ are required. If all the saddle points leading to different system states can be located, the activation energies of each of these pathways can be supplied to the KMC procedure and, with the pre-Arrhenius frequency factor, the rates of each pathway determined. The system can then be propagated in a dynamically correct way to the next state, and the entire procedure is repeated.

In practice, however, it is computationally demanding (effectively impossible) to find all the saddle points and demonstrate that all have been found, and the number of saddle points bounding a state grows exponentially with the dimensionality of the system. Xu and Henkelman²² have developed a method to quantify a ‘confidence limit’ for each state that enough relevant saddle points have been found to progress the simulation with a pre-defined level of confidence and this approach has been implemented here. Their method supposes that all of the relevant kinetic events possess activation energies within the range mk_BT of the lowest barrier process. m should be carefully chosen, and large enough that relevant transitions are not omitted. For example, for m = 20, an event at the upper limit will be e⁻²⁰ ≈ 10⁻⁹ as likely to occur as one at the lower limit, assuming that the pre-factors are equal. Events that occur within these limits are therefore relevant to the kinetics of the system, and the confidence (C) that a relevant saddle has not been missed is:

N_r is the number of sequential searches that find relevant, but redundant (non-unique) processes and α lies between zero and one and describes the relative probability of finding each relevant saddle point. α = 1 describes a system where there is an equal probability of finding each relevant saddle point, and is the value used in this work. More justification of eqn (1) is detailed elsewhere.²² By setting a confidence limit of, e.g. 95%, the simulation will run until N_r = 20 searches finish without finding a new, unique saddle point. Each saddle point search proceeds assuming no knowledge of the local environment, with no saddle point recycling.

We use a new adaptive KMC program, DL_AKMC,† as follows. Initially the provided material structure is minimised using a user-specified algorithm as implemented in DL_FIND.²³ Possibilities are steepest descent, conjugate gradient following Polak–Ribiére,²⁴ L-BFGS,^25,26 P-RFO,^27–30 Newton–Raphson/quasi-Newton, damped dynamics, random (stochastic) search^31,32 or by a genetic algorithm.^33–35 This minimised structure is then used as the initial basin for the task-farmed saddle-point searches either by the improved dimer method^19,20 or the NEB method³⁶ as implemented in DL_FIND, as specified by the user. These dimer searches are initiated near the local minimum by displacing the system away from the minimum in a random direction. The algorithm for displacement used here is general – each atom is displaced up to a maximum of ±0.4 Å.

The kinetic Monte Carlo algorithm has been documented elsewhere^37,38 and Fig. 1 gives a brief overview of the parallelisation used. Once all of the transitions have been identified, the rate ν_i for each event i is then obtained from eqn (2):

ν₀ is the pre-exponential factor, here set equal to a typical vibrational frequency (10¹³ s⁻¹) as is usual in the absence of more detailed information.¹² ΔE is the activation energy barrier for ion migration determined using the dimer method, k_B the Boltzmann constant and T the temperature. One transition is chosen at random with a probability proportional to the relative rate of the transition. The timestep of the simulation is then advanced by Δt (eqn (3)):ω is a random number between zero and one. The system structure is then updated to reflect the transition chosen, and the saddle-point searching algorithm begins again. This process continues until the maximum simulation or computational time is reached.

Fig. 1 Schematic of the parallelisation and simulation process implemented in DL_AKMC.

2.2 Ionic conductivity

Key to the operation of an SOFC electrolyte or cathode is its ionic conductivity. For efficient functioning of the fuel cell both electrolytes and cathodes should possess high such conductivity. This is related to the oxygen self diffusion coefficient by the Nernst–Einstein relation (eqn (4)):σ_i is the ionic contribution to the conductivity from species i, z_ie the charge of species i, c_i the concentration of ionic defects (in this paper, oxide-ion vacancies) and D_i the diffusion coefficient of species i. This approach is strictly only valid for dilute systems, and not our doped systems which are almost always concentrated solid solutions.

Grope et al.³⁹ and Pornprasertsuk et al.¹¹ have simulated such systems circumventing this issue by direct inclusion of an applied field. In their approach, the electrolyte is split into a series of ‘slices’ that are infinite in two dimensions (y- and z-) perpendicular to the direction of the applied field (see Fig. 2). The applied field alters the activation energy barrier for ion migration (eqn (5)):

ΔE^corr = ΔE⁰ + αz_ieV_shift (5)

ΔE^corr is the new, corrected activation energy barrier for ion migration, ΔE⁰ the activation energy barrier obtained from the saddle-point searches, V_shift is the total potential difference between the destination slice and the initial slice in which the vacancy is located. α is termed the ‘symmetry factor’ and has a value between zero and one which is dependent on the profile of the migration barrier under the applied potential but in this work is set to 0.5. For any oxide-ion vacancy movement that occurs within a slice, V_shift is zero. The total potential, Vⁱ_total, at each slice, i, is split into two terms – the electrode potential (Vⁱ_elec) and the space charge potential (Vⁱ_sc):

Vⁱ_total = Vⁱ_elec + Vⁱ_sc (6)

Fig. 2 Generalised illustration of the supercell used in this work, where the ‘slices’ are infinite in y- and z-dimensions and the electric field is applied along the x-direction.

V ⁱ _elec is the potential on slice i due to charge accumulation at the electrodes. Electrons are assumed to have a much higher mobility than oxide ions and so Vⁱ_elec is equal to the difference between Vⁱ_total and Vⁱ_sc. In our simulations, V_elec is assumed to vary linearly across the electrolyte:

V ⁰_elec, the electrode potential in the initial slice, labelled 0, is set to zero for convenience. V^N_elec is the electrode potential on slice N, the final slice in our configuration. The space charge potential per slice, i is given by:

V ⁰_sc, the space charge potential in the initial slice, 0, set to zero. The slice spacing of all the cubic unit cells of the crystal structures considered in this work is half the lattice parameter, a. In using eqn (8) during the simulation, ionic movement is assumed not to change the geometry of the simulation cell. E_i is the electric field associated with each infinite slice, i given by:

ε₀ is the vacuum permittivity and ε_r the relative permittivity of the material, set at 40 in our simulations, equal to the experimental value for 8 mol% YSZ.⁴⁰ Small variations in the permittivity that occur when simulating different materials do not greatly influence the results. The charge density per slice (ρ_i) is the total defect charge in the slice divided by the surface area of the periodic slice. To simulate an impedance measurement an alternating potential, V_appl, is applied:

V_appl(t) = V₀ cos(2πft) (10)

t is the time, f the frequency and V₀ the amplitude of the applied potential. V^N_total (eqn (6)) is set to −V_appl. Under the applied field the activation energies and subsequently the rates (ΔE and ν_i respectively in eqn (2)) for migration of charged species in the field direction change, and the resulting displacement can be used to calculate the contribution to conductivity from species i, σ_i (eqn (11)):〈X〉 is the mean displacement of species i in the applied field, E_field, c_X the concentration of i and t is the simulation time. We assume no change in charge as the charged species migrate.

For the aKMC calculations the materials were simulated for a minimum of 1 ms, and a summary of the parameters required is given in Table 1.

Table 1 Parameters used in the aKMC simulations

Vacuum permittivity, ε0/(F m^−1):	8.854 × 10^−12
Relative permittivity, εr:	40
Vibrational frequency, ν0/Hz:	1 × 10^13
V 0 (eqn (10)/V:	0.5
Frequency of applied field, f/Hz:	1 × 10^7

---

2.3 Application to oxides

We have investigated four solid oxide fuel cell (SOFC) electrolyte materials with the fluorite structure: yttria- and calcia-stabilised zirconica (YSZ and CSZ respectively, with dopant concentrations varying from 2–18 mol%) and gadolinium- and samarium doped ceria (GDC and SDC respectively, with dopant concentrations varying from 5–25 mol%), and one cathode material with the perovskite structure, strontium-doped LaCoO₃ (LSCO, with dopant concentrations varying from 5–80 mol%). Pictures of the two structures are shown in Fig. 3. YSZ, CSZ, SDC and GDC are all stable at intermediate temperature SOFC operating temperatures (up to 900 °C), and YSZ, SDC and GDC all possesses high ionic conductivities (∼10⁻² S cm⁻¹ at 700 °C).^41,42 Doping of zirconia involves substitution of some some Zr⁴⁺ ions with Y³⁺ or Ca²⁺ and of ceria replacement of some Ce⁴⁺ with Gd³⁺ or Sm³⁺. To maintain charge neutrality, one oxygen vacancy is created for every two Zr or Ce ions replaced with Y, Gd or Sm ions. When Zr ions are replaced with Ca ions, one oxygen vacancy is created per ion replaced. This gives rise to a considerable number of oxide ion vacancies which greatly increase the ionic conductivity of these materials relative to the undoped systems.

Fig. 3 The (a) fluorite and (b) perovskite crystal structures. Cation (left) and anion (right) sublattices are separated for clarity. The perovskite structure contains A-site cations (large, grey, central atom in figure), and B-site cations (smaller, blue atoms in figure).

2.4 Potential model

For the adaptive kinetic Monte Carlo simulations a rigid-ion model, including 2-body Buckingham potentials is used (eqn (12)) to model the atomic interactions. The potential parameter set is listed in Table 2. A non-Coulombic potential cutoff of 12 Å was used throughout. Each ion is assigned an integer charge, i.e. +2 for Ca and Sr, +3 for Gd, Sm, La and Co, +4 for Zr and Ce and −2 for O.A, ρ and C are constants, r_ij is the distance between ions of type i and j, and z_ie and z_je are the charges of species i and j respectively.

Table 2 Buckingham potential parameters used in the adaptive KMC calculations. Parameters A, ρ, C are defined in eqn 12

Interaction	A/eV	ρ/Å	C/(eV Å^6)	Reference
Ca^2+–O^2−	1090.40	0.3440	0.00	Ref. 43
Sr^2+–O^2−	959.10	0.3721	0.00	Ref. 44
Y^3+–O^2−	1325.60	0.3461	0.00	Ref. 45
Gd^3+–O^2−	1962.74	0.3250	0.00	Ref. 46
Sm^3+–O^2−	1944.44	0.3414	21.49	Ref. 47
La^3+–O^2−	1545.21	0.3590	0.00	Ref. 48
Co^3+–O^2−	1329.82	0.3087	0.00	Ref. 48
Zr^4+–O^2−	1024.60	0.3760	0.00	Ref. 45
Ce^4+–O^2−	1809.68	0.3547	20.40	Ref. 49
O^2−–O^2−	17 428.92	0.1490	27.89	Ref. 46

---

Since our potential set uses a number of potentials from different sources, we have validated the set by checking good agreement between simulated and experimental values for bulk properties (lattice parameter, elastic constants and bulk modulus) of a number of binary oxides using the GULP code.⁵⁰ The results are shown in the ESI, Section 3.†

Cubic supercells containing 4 × 4 × 4 cubic unit cells (for SDC, GDC, YSZ and CSZ) and 6 × 6 × 6 cubic unit cells (for LSCO) were constructed and both the dopant cations and required oxygen vacancies were placed randomly in the lattice, subject to the constraint that each slice (0 to N, see Fig. 2) has zero net charge. It has been suggested that the oxygen vacancies in YSZ preferentially occupy sites near to yttrium ions (second nearest neighbour positions⁵¹), which is not reflected in our choice of initial structure (randomly distributed stabilisers and vacancies). We have performed studies investigating the influence of the starting structure on the ionic conductivities at the end of the adaptive kinetic Monte Carlo simulations. Both structures with anion vacancies and dopant cations as nearest neighbours, and structures with anion vacancies and dopant cations separated by a minimum of 4.5 Å result in ionic conductivities similar to those seen with the random arrangement used here. The initial arrangement of the dopant atoms and vacancies appears to bear little influence on the ionic conductivity of the material at the end of the simulations. More detail on this is given in the ESI, Section 1.† The initial random structures used in this work were equilibrated for 20 ps (NPT ensemble) using DL_POLY⁵² at the temperature used in the subsequent DL_AKMC simulation (typically 300 K, with the exception of those used to create the Arrhenius plot in Fig. 9 later). After equilibration each slice was checked to ensure charge neutrality was maintained, with the resulting arrangements used as input structures for DL_AKMC.

To check further the quality of the potentials employed, the lattice parameters of the simulation cells (after equilibration) are compared with experimental lattice parameters as a function of the dopant concentration (Fig. 4 and 5).

Fig. 4 Calculated (solid lines) and experimental (dashed lines) lattice parameters for YSZ⁵³ (circles, black), CSZ (squares, red), GDC⁴² (triangles, green) and SDC⁴² (diamonds, blue) as a function of dopant concentration. Lines are intended as a guide for the eye.

Fig. 5 Calculated (circles, black solid line) and experimental⁵⁴ (squares, red dashed line) lattice parameters for LSCO as a function of strontium content. Lines are intended as a guide for the eye.

Fig. 4 shows that the calculated lattice parameters for the doped zirconias and cerias agree with the experimental values with a match of <1%. Similarly, the lattice parameters for LSCO are in good agreement, with the calculated lattice parameter slightly larger then experiment (by 0.6–1.6%). Both experiment and simulation agree as to the increase of the lattice parameter of the doped system with increasing dopant content. Our potentials thus describe these systems adequately and are suitable for use with adaptive kinetic Monte Carlo.

3 Results and discussion

Conductivities were calculated (eqn (11)) for all the materials and the results are shown in Fig. 6 and 7 for the electrolyte materials and Fig. 8 for LSCO, along with some representative experimental values. It can be seen that there is good agreement between the simulations and experiment for all of the materials modelled. In particular the calculated values of the ionic conductivities in the doped cerias (Fig. 7) match very closely the experimental magnitude of the conductivity, and its variation with dopant concentration. The calculated maximum in conductivity for SDC occurs at 16 mol%, higher than that found experimentally by 2 mol%. Our choice of electrolytes has been governed by their high ionic conductivity and low electronic conductivity and so we would expect that our simulations, which solely determine the ionic conductivity, should agree with the experimental values which are the sum of ionic and electronic contributions. The good agreement found between the calculated and experimental conductivities is emphasised by our use of linear y-axes, in contrast to the logarithmic scale typically deployed for such figures.

Fig. 6 Ionic conductivities for YSZ (calculated: black solid line, circles; experimental:⁶⁰ black dashed line, squares) and CSZ (calculated: red solid line, diamonds; experimental:⁵⁸ red dashed line). Lines are intended as a guide for the eye.

Fig. 7 Ionic conductivities for SDC (calculated: black solid line, circles; experimental:⁴² black dashed line, triangles) and GDC (calculated: red solid line, squares; experimental:⁴² red dashed line, diamonds). Lines are intended as a guide for the eye.

Fig. 8 Calculated (solid black line, circles) and experimental (red squares,⁵⁹ diamonds⁶¹ and triangles⁶²) ionic conductivities for LSCO. The line is intended as a guide for the eye.

The form of the ionic conductivity versus dopant concentration figures has been discussed by others.^55,56 For YSZ in particular, the decrease in conductivity at larger dopant concentrations has been attributed, at least in part, to the formation of strongly-bound complexes reducing the availability of free available for diffusion.⁵⁷ CSZ exhibits a more unusual curve, with two maxima appearing in the dopant range simulated. While the location and magnitude of the second maximum at 14 mol% matches well with the experimental data available,⁵⁸ the general paucity of experimental ionic conductivity measurements for CSZ is such that there are no data for comparison at 6 mol%, and so the first maximum remains an interesting prediction for future verification. It is possible this first maximum may be simply an artifact of the simulation, and a discussion on the reproducibility of results and influence of the starting configuration is given in the ESI, Sections 1 and 2.†

For a mixed electronic and ionic conductor such as LSCO the electronic conductivity σ_el ∼ 10²–10³ S cm⁻¹, dominates.⁵⁹ The simulated values of the ionic conductivity (σ_ion) shown in Fig. 8 cannot be compared directly with experiment. However values of σ_ion have been calculated from experimental oxygen diffusion coefficients and the Nernst–Einstein equation. These values are in the range 10⁻³ to 10⁻¹ S cm⁻¹ at 900 °C depending on composition.⁵⁹ These experimentally-derived values, relate to different temperatures (1048–1173 K), but do highlight the broad margin of error inherent in these measurements. Our theoretically determined values fall within the experimental range of ionic conductivities, and exhibit a coherent trend across the entire dopant range. A maximum in the ionic conductivity is found at ∼48 mol% Sr. Combined with the good agreement between theory and experiment for the electrolytes, we are confident that our values for the ionic conductivity of the LSCO cathode are accurate. Materials such as LSCO, where experimental values are difficult to obtain, emphasise the importance of these simulations for extracting important physical properties that would otherwise be difficult or impossible to measure experimentally.

For each material, one experimentally common dopant concentration was chosen (10 mol% for YSZ and CSZ, 20 mol% for SDC and GDC and 50 mol% for LSCO) and simulated over a temperature range of 600–1100 K. The conductivities were determined over these temperature ranges and the resultant Arrhenius plot is shown in Fig. 9. The average activation energy of the migration processes occurring in each material can be determined from the gradient of a linear fit to the data. For the stabilised zirconias this average activation energy is 0.46 ± 0.03 eV for YSZ and 0.44 ± 0.03 eV for CSZ. The YSZ activation energy, while lower than the experimentally determined values^63,64 of 0.8–1.0 eV, is in keeping with other simulated values^65,66 which range from 0.2–0.9 eV. The doped cerias exhibit a lower average activation energy for diffusion, with calculated values of 0.23 ± 0.04 eV for SDC and 0.28 ± 0.03 eV for GDC, lower than the experimentally determined values of around 0.9 eV for 20 mol% SDC and GDC.⁶⁷ LSCO has an average activation energy similar to that of the doped cerias, 0.32 ± 0.02 eV. Experimentally for LSCO, reported activation energies of the migration processes^59,68 vary greatly from 0.6–2.2 eV for 20 mol% Sr. De Souza et al.⁵⁹ also found that the activation energy dropped substantially on increasing the strontium content, with the activation energy dropping to 1.4 eV for 50 mol% Sr from 2.2 eV for 20 mol%. In general we find that the calculated activation energies for the migration processes in these materials are lower than those found experimentally. This can possibly be attributed to the use of a fixed charge, rigid-ion model where the ionic charge does not change throughout the simulation; in reality the charges on the ions are likely to change depending on the local environment.

Fig. 9 Arrhenius plots for ionic conductivities of 10 mol% YSZ (black, circles), 10 mol% CSZ (red, squares), 20 mol% GDC (green, diamonds), 20 mol% SDC (blue, triangles) and 50 mol% LSCO (orange, crosses).

4 Conclusion

The reliability of the adaptive kinetic Monte Carlo program DL_AKMC has been verified by the accurate simulation of ionic conductivities of the solid electrolytes yttria- and calcia-stabilised zirconica (YSZ and CSZ respectively), gadolinium- and samarium doped ceria (GDC and SDC respectively), and of one cathode material with the perovskite structure, strontium-doped LaCoO₃ (LSCO). The calculated ionic conductivities of all the materials are of the same order of magnitude as those found experimentally. The variation of ionic conductivity with dopant concentration follows the same form as experiment. The simulations also indicate the presence of a possible peak in the conductivity of CSZ at ∼6 mol%, which would be of interest to probe experimentally. Average activation energies of all migration processes are between 0.4–0.5 eV for the stabilised zirconias, between 0.2 and 0.3 eV for the doped cerias and 0.3 eV for LSCO.

This work has showcased the ability of the program to simulate experimentally measurable and important properties, while retaining the flexibility of using common potential types and the ability to determine transition states on-the-fly in a massively parallel fashion. Measuring solely the ionic conductivity in mixed ionic-electronic conductors is difficult, and there is a large variation in the experimental ionic conductivity, particularly with LSCO. The ability for a simulation to model reliably the ionic conductivity of such mixed conductors with high precision and accuracy is of critical importance, and such predictions are now achievable and efficient. Given the flexibility and general applicability of the DL_AKMC program, it will likely prove to be a vital tool in a wide range of scientific fields.

The authors would like to thank John Harding, Karl Travis and Ilian Todorov for helpful conversations. This work has been supported by the HECToR Distributed Computational Science and Engineering (CSE) Service operated by NAG Ltd., and by the UK Engineering and Physical Sciences Research Council (EPSRC) under Grant EP/H012990/1.

References

1. T. H. Etsell and S. N. Flengas, J. Electrochem. Soc., 1890, 118, 1971.
2. D. P. Fagg, V. V. Kharton and J. R. Frade, J. Electroceram., 2002, 9, 199.
3. T. S. Zhang, J. Ma, Y. J. Leng, S. H. Chan, P. Hing and J. A. Kilner, Solid State Ionics, 2004, 168, 187.
4. F. M. Figueiredo, F. M. B. Marques and J. R. Frade, Solid State Ionics, 1998, 111, 273.
5. A. Y. Zuev, A. N. Petrov, A. I. Vylkov and D. S. Tsvetkov, J. Mater. Sci., 1901, 42, 2007.
6. W. Zhou, W. Jin, Z. Zhu and Z. Shao, Int. J. Hydrogen Energy, 2010, 35, 1356.
7. H. Fukunaga, M. Koyama, N. Takahashi, C. Wen and K. Yamada, Solid State Ionics, 2000, 132, 279.
8. Z. P. Shao and S. M. Haile, Nature, 2004, 431, 170.
9. R. Sayers, N. L. O. Flack, J. Alaria, P. A. Chater, R. G. Palgrave, S. R. C. McMitchell, S. Romani, Q. M. Ramasse, T. J. Pennycook and M. J. Rosseinsky, Chem. Sci., 2013, 4, 2403.
10. P. F. Shimojo and H. Okazaki, J. Phys. Soc. Jpn., 1992, 61, 4106.
11. R. Pornprasertsuk, J. Cheng, H. Huang and F. B. Prinz, Solid State Ionics, 2007, 178, 195.
12. R. Krishnamurthy, Y. G. Yoon, D. J. Srolovitz and R. Car, J. Am. Ceram. Soc., 1821, 87, 2004.
13. D. J. Harris, T. S. Farrow, J. H. Harding, M. Y. Lavrentiev, N. L. Allan, W. Smith and J. A. Purton, Phys. Chem. Chem. Phys., 1839, 7, 2005.
14. M. Y. Lavrentiev, D. J. Harris, J. H. Harding, N. L. Allan and J. A. Purton, Dalton Trans., 2004, 3071.
15. D. J. Harris, M. Y. Lavrentiev, J. H. Harding, N. L. Allan and J. A. Purton, J. Phys.: Condens. Matter, 2004, 16, L187.
16. G. Henkelman and H. Jónsson, J. Chem. Phys., 2001, 115, 9657.
17. C. Wert and C. Zener, Phys. Rev., 1949, 76, 1169.
18. G. H. Vineyard, J. Phys. Chem. Solids, 1957, 3, 121.
19. G. Henkelman and H. Jónsson, J. Chem. Phys., 1999, 111, 7010.
20. J. Kästner and P. Sherwood, J. Chem. Phys., 2008, 128, 014106.
21. A. F. Voter, Phys. Rev. Lett., 1997, 78, 3908.
22. L. Xu and G. Henkelman, J. Chem. Phys., 2008, 129, 114104.
23. J. Kästner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander and P. Sherwood, J. Phys. Chem. A, 2009, 113, 11856.
24. E. Polak and G. Ribiére, Ref. Fr. Inf. Rech. Op., 1969, 3, 35.
25. D. C. Liu and J. Nocedal, Math. Program., 1989, 45, 503.
26. J. Nocedal, Math. Comp., 1980, 35, 773.
27. C. J. Cerjan and W. H. Miller, J. Chem. Phys., 1981, 75, 2800.
28. J. Simons, P. Jørgensen, H. Taylor and J. Ozment, J. Phys. Chem., 1983, 87, 2745.
29. A. Banerjee, N. Adams, J. Simons and R. Shepard, J. Phys. Chem., 1985, 89, 52.
30. J. Baker, J. Comput. Chem., 1986, 7, 385.
31. S. H. Brooks, Oper. Res., 1957, 6, 244.
32. R. Luus and T. H. I. Jaakola, AIChE J., 1973, 19, 760.
33. J. H. Holland, Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, 1975.
34. D. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, New York, 1989.
35. R. L. Haupt and S. E. Haupt, Practical genetic algorithms, Wiley, New York, 1998.
36. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978.
37. A. B. Bortz, M. H. Kalos and J. L. Lebowitz, J. Comp. Physiol., 1975, 17, 10.
38. D. T. Gillespie, J. Comput. Phys., 1976, 22, 403.
39. B. O. H. Grope, T. Zacherle, M. Nakayama and M. Martin, Solid State Ionics, 2012, 225, 476.
40. D. P. Thompson, A. M. Dickins and J. S. Thorp, J. Mater. Sci., 1992, 27, 2267.
41. P. K. Moon and H. L. Tuller, Solid State Ionics, 1988, 28, 470.
42. S. Zha, C. Xia and G. Meng, J. Power Sources, 2003, 115, 44.
43. N. H. de Leeuw, G. W. Watson and S. C. Parker, J. Phys. Chem., 1995, 99, 17219.
44. G. V. Lewis and C. R. A. Catlow, J. Phys. C: Solid State Phys., 1985, 18, 1149.
45. M. Kilo, C. Argirusis, G. Borchardt and R. A. Jackson, Phys. Chem. Chem. Phys., 2003, 5, 2219.
46. D. S. D. Gunn, N. L. Allan, H. Foxhall, J. H. Harding, J. A. Purton, W. Smith, M. J. Stein, I. T. Todorov and K. P. Travis, J. Mater. Chem., 2012, 22, 4675.
47. L. Minervini, R. W. Grimes and K. E. Sickafus, J. Am. Ceram. Soc., 1873, 83, 2000.
48. M. Cherry, M. S. Islam and C. R. A. Catlow, J. Solid State Chem., 1995, 118, 125.
49. S. Vyas, R. W. Grimes, D. H. Gay and A. L. Rohl, J. Chem. Soc., Faraday Trans., 1998, 94, 427.
50. J. Gale, Z. Kristallogr., 2005, 220, 552.
51. H. Ding, A. V. Virkar and F. Liu, Solid State Ionics, 2012, 215, 16.
52. W. Smith, T. R. Forester and I. T. Todorov, The DL_POLY Classic User Manual, Daresbury Laboratory, 2010.
53. C. Pascual and P. Duran, J. Am. Ceram. Soc., 1983, 66, 23.
54. A. N. Petrov, O. F. Kononchuk, A. V. Andreev, V. A. Cherepanov and P. Kofstad, Solid State Ionics, 1995, 80, 189.
55. A. I. Ioffe, D. S. Rutman and S. V. Karpachov, Electrochim. Acta, 1978, 23, 141.
56. S. P. S. Badwal, Solid State Ionics, 1992, 52, 23.
57. R. Devanathan, W. J. Weber, S. C. Singhal and J. D. Gale, Solid State Ionics, 2006, 177, 1251.
58. T. Takahashi, High Conductivity Solid Ionic Conductors, World Scientific, Singapore, 1989.
59. R. A. De Souza and J. A. Kilner, Solid State Ionics, 1998, 106, 175.
60. M. Filal, C. Petot, M. Mokchah, C. Chateau and J. L. Carpentier, Solid State Ionics, 1995, 80, 27.
61. E. Bucher, W. Sitte, I. Rom, I. Papst, W. Grogger and F. Hofer, Solid State Ionics, 2002, 152, 417.
62. Y. Teraoka, T. Nobunaga, K. Okamoto, N. Miura and N. Yamazoe, Solid State Ionics, 1991, 48, 207.
63. R. E. W. Casselton, Phys. Status Solidi A, 1970, 2, 571.
64. J. Luo, D. P. Almond and R. Stevens, J. Am. Ceram. Soc., 2000, 83, 1703.
65. X. Li and B. Hafskjold, J. Phys.: Condens. Matter, 1995, 7, 1255.
66. M. Kilo, R. A. Jackson and G. Borchardt, Philos. Mag., 2003, 11, 3309.
67. V. Esposito and E. Traversa, J. Am. Ceram. Soc., 2008, 91, 1037.
68. S. Carter, A. Selcuk, R. J. Chater, J. Kajda, J. A. Kilner and B. C. H. Steele, Solid State Ionics, 1992, 53, 597.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ta01504e

---