ARTICLE pubs.acs.org/JPCC
Atomistic Simulation of the Structural, Thermodynamic, and Elastic Properties of Li2TiO3 Samuel T. Murphy,*,†,‡ Philippe Zeller,† Alain Chartier,† and Laurent Van Brutzel† † ‡
CEA, DEN, LM2T, F-91191 Gif-sur-Yvette, France Department of Materials, Imperial College London, South Kensington, London, SW7 2AZ, U.K. ABSTRACT: Lithium-based ceramics, such as lithium metatitanate, have been proposed for adoption in the breeder blanket region of a fusion reactor. In this article, we report a combination of empirical and density functional theory (DFT) simulations employing “on-the-fly” pseudopotentials for Li2TiO3. The smoothing parameters of the planewave pseudopotentials were optimized to ensure an appropriate level of precision for determination of structural, thermodynamic, and elastic properties. As the elastic properties of lithium metatitanate are not well-known, the efficacy of the DFT simulations employing the new pseudopotentials was explored using Li2O and TiO2 where experimental data are available. These pseudopotentials are then used to investigate the three intermediate temperature phases of Li2TiO3 (i.e., C2/c, C2/m, and P3112). Finally, we examine the elastic properties of Li2TiO3 using both DFT and an empirical potential model and find it to be, irrespective of space group, more resistant to deformation than other promising ceramic breeder materials.
’ INTRODUCTION Current efforts to realize the energy generating potential of nuclear fusion include the tokamak, in which a fusion plasma is suspended in a toroidal magnetic field. Within the fusion plasma, deuterium and tritium are fused together to form helium, and a highly energetic neutron is ejected (see eq 1). Deuterium can be easily extracted from seawater; however, due to its relatively short halflife, there is no naturally available source of tritium, and it is expedient to generate it in situ from the transmutation of lithium (see eq 2 and eq 3).1 2 1D
þ 31 T f 42 He þ 10 n
ð1Þ
6 3 Li
þ 10 n f 42 He þ 31 T
ð2Þ
7 3 Li
þ 10 n f 42 He þ 31 T þ 10 n
ð3Þ
The production of tritium via the transmutation of lithium will occur in the breeder blanket region of a fusion reactor. Situated immediately behind the plasma facing first wall, the breeder blanket will serve three main purposes: (i) the production of tritium as already mentioned, (ii) conversion of the neutron energy into heat energy that can be used in electricity generation (this is because the majority of the energy released during the DT reaction is localized on the neutron in eq 1),2 and (iii) to protect the magnets and the superstructure from neutron and gamma irradiation. To maintain a sustainable and economic process, it is essential that for every neutron ejected from the plasma a tritium atom can be recovered from the breeder blanket region (i.e., the tritium breeding ratio, TBR, must be greater than one). A high TBR can be achieved by selecting blanket materials r 2011 American Chemical Society
with high lithium densities and enrichment of lithium to ensure the neutron is recycled (eq 3), as well as the introduction of neutron multiplier materials such as beryllium.3 Due to their high melting temperatures, low chemical reactivity, and high lithium densities, lithium-containing ceramic oxides are an attractive group of materials for adoption in the breeder blanket region of a fusion reactor. Candidate ceramic breeder materials include lithium metatitanate (Li2TiO3),4 lithium orthosilicate (Li4SiO4),5 and lithium oxide (Li2O).4 In addition to possible application in fusion reactors, the high lithium mobility demonstrated by these materials has seen them applied as electrode materials in lithium-ion batteries.6 There has been significant interest over the past decade in generating a complete database of relevant properties of these lithium ceramic materials. Despite this, for lithium metatitanate, there are still gaps in the database; some mechanical properties remain only partially known (in particular, the elastic moduli are either unknown or the data display a wide distribution7,8), and tritium migration through the matrix is also poorly understood. The recent growth in computational power has enabled the simulation of large periodic systems using ab initio methods, such as density functional theory (DFT). One of the major approximations adopted to ensure DFT calculations are computationally tractable is the introduction of pseudopotentials. Pseudopotential theory assumes that the quantum state of core electrons is invariant under the various chemical bonding and/or physical conditions, such as pressure, to which the element will be Received: May 19, 2011 Revised: September 20, 2011 Published: September 28, 2011 21874
dx.doi.org/10.1021/jp204678c | J. Phys. Chem. C 2011, 115, 21874–21881
The Journal of Physical Chemistry C submitted in further computations. This quantum state is computed in the plane-wave pseudopotential construction phase. The number of electrons encapsulated in the core is a crucial design choice: a smaller core offers a more accurate description of compounds but at a higher computational cost. This trade-off applies more generally to most other pseudopotential design parameters. Pseudopotentials are often taken from a standard library, which may be either linked with the DFT code itself or independent. In both cases these libraries are generally optimized for versatility; i.e., each element is treated separately to give an average performance with all other elements. Here we explore the possible benefits of using the “on-thefly” pseudopotential generating feature in CASTEP that allows customization of the pseudopotential specific to the system investigated and the specific properties to be studied. Many of the atomistic simulation works on lithium ceramics have focused on the binary lithium oxide,913 and there is a surprising paucity of data for the ternary oxides. Zainullina et al.14 used the linear “muffin-tin” orbitals method to investigate the incorporation of hydrogen on lithium sites and predict substitution preferentially on the lithium site in the LiTi2 layer (see eq 3). HartreeFock cluster simulations15,16 have suggested that tritium would preferentially occupy the octahedral interstitial site close to the pure lithium layer as opposed to substitution onto a lithium site. The diffusion, and ultimately the release rate, of tritium may be linked to the diffusion on the lithium sublattice. To investigate this, Vijayakumar et al.17 performed molecular dynamics simulations (coupled with 6,7 Li NMR studies) to investigate the diffusion of Li+. Their results suggest that a lithium vacancy will preferentially reside in the LiTi2 layer and that the effective activation energy for vacancy migration is 0.30 eV. Overall, they found that lithium ion conductivity for Li2TiO3 is less than for other traditional lithium ion conductors. Density functional theory has also been used to study the relative energies of lithiated TiO2 polymorphs.18 The objectives of this article are 3-fold: (i) to optimize the smoothing parameter for the future study of systems such as Li2TiO3, (ii) to demonstrate the efficacy of DFT employing the new pseudopotentials by providing a rigorous validation for Li2O and TiO2 where more experimental data are available, and (iii) to use this model in addition to an empirical model developed by Vijayakumar et al.17 to investigate the complex structural, thermodynamic, and elastic properties of lithium metatitanate.
’ CRYSTALLOGRAPHY The Li2OTiO2 phase diagram19,20 shows that Li2TiO3 crystallizes in three phases. The metastable cubic α-phase transforms to the monoclinic β-phase at 300 °C. This is the stable phase until 1215 °C after which it transforms to the cubic γphase. A further phase transition was predicted to occur at 450 °C,21 although the nature of this modification is not described. The crystal structure of β-phase Li2TiO3 was first determined by Lang22 and subsequently refined using X-ray diffraction of large single crystals by Kataoka et al.23 These works envisaged Li2TiO3 as a distorted rocksalt structure with alternating Li and LiTi2 (111) planes described with the C2/c space group. More recently, Tarakina et al.24 examined the structure of Li2TiO3 at 700 °C and found that the intensity of some Bragg reflections and broad asymmetric diffuse scattering in the 2θ range 20°30° of their electron diffraction patterns could not be refined using an ideal model based on the C2/c space group. They proposed a model of defect stacking of the LiTi2 layers similar to that observed in Li2MnO3.25 The arrangement of the
ARTICLE
Figure 1. LiTi2 layer stacking schemes in Li2TiO3. (a) Shows the honeycomb LiTi2, A1, layer where green spheres represent the lithium ions, yellow spheres represent titanium atoms, and black spheres represent the possible stacking locations for the central lithium ions. (b), (c), and (d) show the different stacking schemes in the z-direction that generate the C2/m, P3112, and C2/c space groups, respectively.
atoms in each LiTi2 layer is identical and can be described as a honeycomb structure in the xy plane with six Ti atoms creating a hexagon centered on a Li atom. Figure 1(a) shows an example of this honeycomb structure. The subsequent stacking sequences of these honeycomb structures then give rise to a number of different space groups for Li2TiO3. When a second honeycomb layer is stacked onto the first, there is only one stacking possibility, A1B1 (all other arrangements that maintain the close-packed structure will be related by rotation over 120°). There are two possibilities for addition of a third honeycomb layer: A1B1C1 (Figure 1(b)), which introduces a mirror plane resulting in the C2/m space group, and A1B1C2 (Figure 1(c)), which introduces a 3-fold screw axis along the zdirection giving the P3112 space group. More complex stacking regimes are possible if these three layers are combined with others. For example, combining A1B1C2 with A2B3C1 results in the six-layer A1B1C2A2B3C1 illustrated in Figure 1(d). This stacking arrangement results in a glide plane and the C2/c space group as predicted by Kataoka et al.23 Using the crystal structures discussed above, Tarakina et al.24 controlled the probability of each space group to appear during the generation of stacked layers and prepared a theoretical X-ray diffraction pattern. This theoretical X-ray diffraction pattern was then compared to their 700° sample. After refinement, their model yielded probabilities for short-range ordering of 34.7%, 46.6%, and 19.2% for C2/c, C2/m, and P3112 structures, respectively.
’ COMPUTATIONAL PROCEDURE Density Functional Theory. All DFT simulations presented here were conducted using the CASTEP 5.0 simulation package.26,27 CASTEP is a plane-wave pseudopotential code that describes a crystal using supercells and periodic boundary conditions with special point 21875
dx.doi.org/10.1021/jp204678c |J. Phys. Chem. C 2011, 115, 21874–21881
The Journal of Physical Chemistry C
ARTICLE
Figure 2. Plots showing the difference in energy (a) and the log of the difference in energy (b) between Li, Ti, and O determined at a given cutoff energy and the same atom estimated for an infinite cutoff energy as a function of the plane-wave cutoff energy.
rij 1 qi qj þ Aij exp Uij ¼ 4πε0 rij Fij
!
Cij rij6
ð4Þ
) )
) )
where ε0 is the permittivity of free space; Aij, Fij, and Cij are shortrange potential parameters specific to ions i and j; rij is the interionic separation; and qi and qj are the charges on i and j, respectively. Here we employ the potential model of Vijayakumar et al.17 whereby the ions are assigned noninteger, rather than formal, charges (i.e., qLi = 0.549 e , qTi = 2.196 e , and qO = 1.098 e ). The potential parameters were derived in stages with the TiO parameters originally developed by Matsui and Akaogi32 (MA) by fitting to the lattice parameters of the TiO2 polymorphs: rutile, anatase, Brookite, and TiO2II. Due to the large diffuse nature of the oxygen ions, Kerisit et al.33 added a shell model, such that a massive positively charged core is attached to a negatively charged shell by means of a harmonic force constant, to the original MA potentials by fitting to DFT simulations. Finally the LiO parameters were developed by fitting to the experimental lattice parameters and bulk modulus of Li2O.17 Importantly, these potential parameters were not derived specifically for the Li2TiO3 system. All empirical simulations conducted in the study were performed using the GULP simulation package.34 ) )
integration over the Brillouin zone. We employ the generalized gradient approximation of Perdew, Burke, and Ernzerhof (GGA PBE).28 A MonkhorstPack29 scheme is used to sample the Brillouin zone with a density of 0.05 Å1 on each axis. The Fourier transform grid for the electron density is larger than that of the wave functions by a scaling factor of 2.0 (represented in CASTEP using the grid_scale parameter), and the corresponding scaling for the augmentation values (the fine_grid_scale parameter) was set to 2.3. These values were determined by performing convergence tests of the energy from self-consistent single-point simulations. To ensure that the optimum precision was obtained, we use a convergence criterion for the self-consistent simulations of 5 106 eV. Similarly, robust criteria for the geometry optimization (energy = 5 105 eV atom1, forces = 0.01 eV/Å, stress = 0.02 GPa, and displacement = 5 104 Å) were employed. In this study, we make use of the “on-the-fly” pseudopotential (PP) generating feature to generate a series of ultrasoft pseudopotentials, based on the formalism of Vanderbilt,30 as implemented in CASTEP. Upon initialization CASTEP performs a reference, all electron (AE) calculation on a user-defined atomic configuration at the relevant level of theory (as mentioned previously here we use GGA-PBE28). Pseudo wave functions (PWF) are then fitted to the all electron valence orbitals (AEVO) outside a defined cutoff radius, rc. The PP is made such that the PWF and augmentation charges are as smooth as possible within rc while mimicking the AEVO outside of the core region and ensuring the eigenenergies are as close as possible to the AE eigenenergies. There is, however, a trade-off between the pseudodensity smoothness and the closeness of the PWF and pseudoenergies to all electron ones. Smoothing of the PWF is controlled in CASTEP using the qc optimization scheme (i.e., a low qc results in a soft PP, while a higher value leads to a harder PP with an improved level of agreement between the PWF and the AE wave functions). It is this parameter we seek to optimize here. Empirical Simulations. The empirical simulations conducted in this study envisage the crystal as an infinite array of point charges as first described by Born.31 Interactions between these point charges are then represented using a combination of a long-range Coulombic interaction and a short-range isotropic Buckingham potential (containing a van der Waals contribution), as shown in eq 4
’ RESULTS AND DISCUSSION Pseudopotentials. The pseudopotentials were generated by assuming the atoms are in their ground states. Semicore 3s and 3p electrons were included as valence when generating the pseudopotential for Ti. We arbitrarily set the desired level of precision between the energy of a single atom determined using the new PP and the energy that would be obtained with an infinite plane-wave cutoff energy to Et < 102 eV. The starting points for the procedure are the default “on-the-fly” parameters in CASTEP. For each atom (i.e., Li, Ti, and O), the parameter controlling the hardness (or smoothness), qc, of the pseudopotential was varied to minimize the cutoff energy at which Et is obtained. The minimum cutoff energy at which Et < 102 eV could be satisfied for oxygen ≈550 eV, making this the hardest of the three pseudopotentials. Given that a single plane-wave cutoff energy is employed in a simulation, applying to all atoms, we set this cutoff to 550 eV. The qc values for Li and Ti were increased, thus generating harder pseudopotentials, to maximize the precision obtained for this cutoff energy. The results are shown in Figure 2, where the difference in energy between a single atom computed using the new pseudopotential and an estimate of the energy for an infinite cutoff energy is plotted as a function of the cutoff energy. Besides illustrating that the O pseudopotential is 21876
dx.doi.org/10.1021/jp204678c |J. Phys. Chem. C 2011, 115, 21874–21881
The Journal of Physical Chemistry C
ARTICLE
Table 1. Comparison of the Lattice Parameters, Elastic Constants, and Bulk Modulus, K, of Li2O and Rutile TiO2 Determined Using DFT and the Customized Pseudopotential and the Empirical Model of Vijayakumar et al.17 (Emp) Compared with Experimental Dataa structure Li2O
parameter a/Å
TiO2
4.62
Emp
experiment
4.54b
4.5435
c11/GPa c12/GPa
203.4 16.5
130.4 67.8
217.035 22.635
c44/GPa
56.2
45.3
68.135
K/GPa
Figure 3. Plot showing the convergence of the total energy as a function of the plane-wave cutoff energy for our pseudopotential and the standard library pseudopotentials. ΔE is taken as the difference between the single-point energy of a nonequilibrium Li2TiO3 with the C2/c space group at cutoff energy, CE, indicated on the x-axis and the single-point energy determined at 650 eV. This plot shows that using our pseudopotential we have obtained an order of magnitude greater precision at a CE of 550 eV than can be obtained using the standard library pseudopotentials.
DFT
a/Å b/Å
78.8 4.64 2.96
b
88.035
88.7
b
4.6036
b
2.9636
4.50 3.00
c11/GPa
276.3
316.2
280.737
c12/GPa
148.2
226.7
189.537
c23/GPa c33/GPa
143.5 468.1
151.7 438.3
155.037 510.837
c44/GPa
122.2
117.5
129.437
c66/GPa
211.6
222.2
211.237
K/GPa
205.3
236.5
210.038
a
the hardest of the three, Figure 2 shows that the precision obtained for Li and Ti is of the order of 103 eV. To demonstrate the convergence achieved using the customized pseudopotentials we calculate the single-point energy per formula unit for a nonequilibrium Li2TiO3 material with the C2/ c space group as a function of the cutoff energy. Figure 3 shows a plot of log[ΔE] as a function of cutoff energy (CE), where ΔE = SP SP ESP CE E650, where ECE is the single-point energy for a given CE SP and E650 is the single-point energy at 650 eV determined using the same PSP. Figure 3 shows that the library pseudopotentials only provided a coarse level of convergence in this cutoff energy range (
