Initial Assessment of an Empirical Potential as a Portable Tool for

Apr 30, 2014 - School of Physics, University of Sydney, Sydney, New South Wales 2006, Australia. J. Phys. Chem. C , 2014, 118 ... Phone: +61 (0)2 9717...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Initial Assessment of an Empirical Potential as a Portable Tool for Rapid Investigation of Li+ Diffusion in Li+‑Battery Cathode Materials Ramzi Kutteh*,†,‡ and Maxim Avdeev† †

Bragg Institute, Australian Nuclear Science and Technology Organization, Lucas Heights, New South Wales 2234, Australia School of Physics, University of Sydney, Sydney, New South Wales 2006, Australia



ABSTRACT: Substantial research activity is currently invested in the pursuit of next generation cathode materials for rechargeable Li-ion batteries. We carry out an initial assessment of the suitability of a recently described empirical potential [J. Phys. Chem. B 2006, 110, 11780] as a rapid, portable, and at least qualitatively accurate computational tool for screening large numbers of potential cathode materials for favorable Li-ion transport capabilities. Selected materials can then be examined more elaborately with more accurate but computationally more expensive first-principles approaches. As test systems for our initial assessment, we chose the group of phosphate olivines LiMPO4 (M = Mn, Fe, Co, Ni), promising candidates for next generation cathode materials and subject of numerous experimental and computational studies. To conduct the assessment, we determined the ground state structures of LiMPO4 from geometry optimizations with this empirical potential and with density functional theory (DFT) and computed activation barriers of Li-ion diffusion in LiMPO4 from molecular dynamics simulations based on the empirical potential and from minimum-energy-path DFT calculations. We show that structural results generated by the empirical potential are in good agreement with the DFT and experimental results and that barrier results produced by this potential are in good agreement with the DFT results and often in better agreement than values generated by custom parametrized empirical potentials.

1. INTRODUCTION Rechargeable batteries are set to play a pivotal role1 in a global economy increasingly environment-aware and reliant on renewable energy resources. Since making their commercial debut more than 20 years ago,2,3 Li-ion rechargeable batteries have rapidly grown in popularity to pervade the power tools and consumer electronics markets, powering in the latter ubiquitous small devices such as cell phones, notebook and tablet computers, and digital cameras, to name only a few. They have largely supplanted in these markets competitors such as NiMH and more conventional NiCd batteries. In the past decade, efforts have been made to scale up Li-ion batteries for applications in the auto industry, serving as power sources for all-electric and plug-in hybrid electric vehicles, as well as in the aviation and space industries. Following the initial reporting4 of a rechargeable Li-metal battery, comprising a lithium anode (negative electrode) and a layered TiS2 intercalation cathode (positive electrode), largely unsuccessful early attempts to address the poor cycling efficiency of lithium in the electrolytic solution (due to the formation of dendrites during cycling, which also presented a safety issue) by alloying it5−10 with other elements (e.g., Al, Si, Sb, Sn, etc.) motivated more successful alternative investigations of intercalation compounds for anode materials,11−13 marking the birth (or rather rebirth14) of the present day “rocking-chair” or, more specifically, Li-ion model of the rechargeable lithium battery. In this model, the Li ions are “rocked” back and forth between the anode and cathode © 2014 American Chemical Society

intercalation materials, such that during the charging process Li ions are de-intercalated from the cathode, which acts as a source of Li ions, travel through the electrolyte, and are intercalated into the anode. During discharge, the reverse process occurs as the Li ions de-intercalate from the anode and intercalate into the cathode, generating a current in an external load. Subsequent crucial introductions of an intercalation LiCoO2 cathode15,16 and an intercalation graphite anode17 eventually culminated in the first commercial and to date most popular rechargeable Li-ion battery type. Li-ion batteries owe their increasing popularity to a range of desirable properties, including high energy density, high voltage, long cycle life, and design flexibility. However, rechargeable Li-ion cell technology is nascent14,18 relative to alternative battery technologies (e.g., lead−acid and NiCd), and improvements in battery performance, calendar and cycle life, safety, and cost are necessary to successfully scale up the technology to the aforementioned heavy equipment and machinery applications (e.g., to make all-electric and plug-in hybrid electric vehicles competitive with conventional combustion engine vehicles), as well as to scale it down to keep pace with the relentless miniaturization of electronic devices. Battery performance in particular is directly related to the characteristics of the electrode materials.14 The lack of current, or even Received: January 15, 2014 Revised: April 22, 2014 Published: April 30, 2014 11203

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

potential next generation, cathode materials with charge storage capacities (e.g., 135 mAh/g for LiCoO2) comparable to those of their anode counterparts (e.g., 372 mAh/g for graphite) has been a major obstacle to improving Li-ion battery performance,19,20 particularly for a successful scale up of the technology to electric vehicles, for example. As a result, a global research effort is underway to find more competitive cathode materials. Here we focus on the computational investigation of Li+ transport in cathode materials and consider as a concrete example the phosphate olivines21 LiMPO4 (M = Mn, Fe, Co, Ni), which have garnered a substantial portion of this research effort. Examining Li+ diffusion in cathode materials by computational means is important for at least two reasons. First, the numerical value of the total (chemical) Li+ diffusion coefficient, which affects the rate capability of the electrode, is capped by the lower of the electronic and ionic conductivities,22 due to the charge neutrality constraint. Hence in electrodes with poor charge/discharge rates, knowledge of the relative magnitudes of the electronic and ionic (i.e., Li+ diffusion coefficient) conductivities can help identify the major reasons behind the poor rate capability. Both high electronic and high Li+ conductivities are needed20 for good rate capability and high power density. Second, from a more fundamental viewpoint, gaining insight into the Li+ diffusion mechanisms and pathways, their relative significance, and how they are impacted by various types of defects, which would be difficult to achieve by purely experimental means, can assist in the design of improved electrode materials (e.g., nanostructured active materials and novel electrode architectures). The olivine phosphates, LiMPO4, form one group of cathode materials potentially more competitive than the layered LiCoO2 used nowadays in most rechargeable Li-ion batteries and thus are currently a subject of intense investigations. Within this group, most of the research activity to date has concentrated on LiFePO4, which has already seen some commercial exposure. Consequently, many computational studies of Li+ diffusion in LiFePO4 have recently appeared, including classical static23,24 and dynamic25,26 investigations using system-fitted customized empirical potentials, and ab initio static27,28 and dynamic29 studies using density functional theory (DFT). Fewer studies30,31 have covered Li+ diffusion in the entire LiMPO4 series. Although our ultimate goal for focusing in this work on the computational investigation of Li+ transport in LiMPO4 is essentially the same as that driving all of these computational studies, namely, to unravel the physics of the Li+ diffusion process in cathode materials as discussed above, our immediate aim for doing so is more intermediary in nature and meant to serve as a stepping-stone to this ultimate goal. Specifically, we seek here a particular computational tool to investigate Li+ diffusion in current and future candidates for next generation Li-ion battery cathode materials. We wish to eventually use this tool to rapidly scan large numbers of potential cathode materials for desirable Li+ transport capabilities and identify those deserving of further detailed investigation with more accurate but also computationally more expensive DFT calculations. Hence the intermediary nature of this desired tool, serving as it were as a link between a database of potential cathode materials and comprehensive Li+ diffusion DFT studies of a select number of these materials, our ultimate goal. To meet our demands, such a tool must (a) produce at least qualitatively correct, compared with DFT calculations, Li+

diffusion coefficients and activation barriers in the cathode materials, (b) do so substantially faster than corresponding DFT calculations, and (c) do so without requiring any systemspecific tuning or reparametrization (i.e., must be portable or transferable among systems). Because this tool is ultimately intended to test a large number of potential cathode materials, requirements a, b, and c must, of course, apply across a wide range of structures and compositions. Rapid computational scanning of materials for particular properties is of course an already familiar paradigm from diverse areas such as virtual screening for drug discovery and design32 and, more recently, high-throughput screening for materials discovery and design.33,34 Our objective at hand is obviously similar in spirit to these endeavors, albeit less ambitious in scope of target materials and properties. In this work, we initiate an assessment of the suitability of an empirical potential, recently described by Pedone et al.35 for silica-based glasses, as a candidate for our desired tool. This potential covers elements including35−38 all of the first row transition metals, as well as Si, P, and O, which makes it particularly relevant for the increasingly popular polyanion cathode materials. A comprehensive assessment of this potential necessarily entails verification of requirements a, b, and c over a representative range of cathode materials, as implied above. Here we embark on this task with a focus on the phosphate olivines, LiMPO4, as test cases. It is worth emphasizing that we are not seeking a computational device for sophisticated studies of Li+ transport in cathode materials (e.g., to examine effects of intrinsic or extrinsic defects on Li+ transport), which are more appropriately performed with DFT or at least with custom parametrized empirical potentials. Our goal is rather a fast, portable, and at least qualitatively accurate precursor to parse through large numbers of cathode compounds and select high Li+ conductivity candidates for more sophisticated DFT calculations, as stated above. Because our choice of computational tool to be assessed here is an empirical potential, it affords by construction calculations of the type mentioned in requirement a above which are much faster than corresponding DFT ones, and hence it automatically satisfies requirement b, irrespective of the chosen test material. Since in addition we have applied this tool throughout the present work to LiMPO4 without system-specific reparametrization, it also satisfies requirement c for this class of materials. Therefore, verification of requirement a for the LiMPO4 test cases is all that remains to complete our initial assessment of this empirical potential. To carry out this verification, we have performed classical calculations of activation energies for Li+ diffusion in LiMPO4 based on the empirical potential of Pedone et al., together with corresponding DFT calculations, and compared the results. These calculations are preceded by computations of the ground state structures of LiMPO4 using this empirical potential and using DFT, with a comparison of the structures obtained with the empirical potential to the DFT and experimental structures partly serving as a prerequisite validation of the empirical potential and its parameters before proceeding to the main task of computing the activation energies. The activation energies generated by the empirical potential are shown to be in good agreement with those from the DFT calculations and, despite requirement c being strictly satisfied, often in better agreement with the DFT results than energies generated by some empirical potentials custom parametrized for LiMPO4. 11204

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

(PAW) method.45,46 Although in its ground state LiMPO4 is antiferromagnetic,47 our DFT calculations were performed conveniently with ferromagnetic spin polarization, the difference in magnetic structure not anticipated to have a major influence on the computed activation energies for Li+ diffusion, as confirmed in section 3.2 by comparison of our activation energies for LiFePO4 with those obtained elsewhere27 using antiferromagnetic spin polarization. It is worth noting that a ferromagnetic ground state for LiMPO4 was also assumed by Morgan et al.30 in their calculations of activation energies of Li+ diffusion. Furthermore, a GGA+U approach was not adopted because, aside from theoretical reservations48 about using this approach to compute diffusion barriers in LiMPO4, it is not our aim here to perform a comprehensive DFT study of Li+ diffusion in LiMPO4 per se, but rather, as discussed in section 1, we wish to obtain DFT activation energies sufficiently accurate to allow a fair comparison with corresponding values generated with the empirical potential of eq 1, verification of requirement a for LiMPO4, and a subsequent initial assessment of this empirical potential. Again, sufficient accuracy of our GGA/DFT activation energies is confirmed in section 3.2 by comparison, for the case of LiFePO4, with activation energies produced elsewhere27 with the GGA+U approach. 2.1. Geometry Optimizations. The olivine crystal structure of LiMPO4 is orthorhombic with space group Pnma.47,49,50 As shown in Figure 1, it comprises a slightly

Section 2 covers the methodologies and parameter values for the structural and Li+ activation barrier calculations of LiMPO4, by both the empirical potential and DFT routes. Section 3 presents and discusses the structural and energy barrier results from these two approaches, together with experimental results and those from other studies based on custom-fitted empirical potentials and DFT. An initial assessment of the empirical potential of Pedone et al. addressing requirement a is given based on this foregoing discussion.

2. COMPUTATIONAL DETAILS We have performed two types of calculations on LiMPO4 (M = Mn, Fe, Co, Ni): ground state structure determination and computation of activation energies of Li+ diffusion along specific pathways in the dilute vacancy limit. Each type was carried out both in the classical domain, based on the empirical potential of Pedone et al.,35 and in the ab initio domain, based on DFT. For each domain, we provide first the computational details common to both types, followed by those specific to each. For the classical calculations on LiMPO4, the empirical potential for a pair of atomic species i and j is given by35,37 ⎧ z z e2 ⎪ i j + D [{1 − e−aij(r − rij̅ )}2 − 1] ij ⎪ r ⎪ ⎪ Cij uij(r ) = ⎨+ 12 , i = Li, M, P, O; j = O ⎪ r ⎪ 2 ⎪ zizje , i , j = Li, M, P ⎪ ⎩ r

(1)

where the first terms in the two expressions above represent together the long-range Coulomb interactions between all pairs of ions in the system, the second term in the top expression represents the short-range i−O interaction (i = Li, M, P, O) modeled by a Morse potential from which the Coulomb interaction has not been subtracted as customary, and the last there is an additional repulsive i−O term required to model interactions at high temperature and pressure. The values35 of the partial charges, zi, and parameters Dij, aij, ri̅ j, and Cij in eq 1 are listed in Table 1. As discussed in section 1, no custom Table 1. Potential Parameters of Equation 1 for LiMPO4 (M = Mn, Fe, Co, Ni)a interaction −1.2

Li −O Mn1.2−O−1.2 Fe1.2−O−1.2 Co1.2−O−1.2 Ni1.2−O−1.2 P3.0−O−1.2 O−1.2−O−1.2 0.6

a

Dij (eV)

aij (Å−2)

ri̅ j (Å)

Cij (eV·Å12)

0.001 114 0.029 658 0.078 171 0.012 958 0.029 356 0.831 326 0.042 395

3.429 506 1.997 543 1.822 638 2.361 272 2.679 137 2.585 833 1.379 316

2.681 360 2.852 075 2.658 163 2.756 282 2.500 754 1.800 790 3.618 701

1.0 3.0 2.0 3.0 3.0 1.0 22.0

Figure 1. Polyhedral depiction of the olivine structure of LiMPO4 (M = Mn, Fe, Co, Ni). The box represents a 1 × 2 × 3 supercell of bulk LiMPO4.

distorted hexagonal close-packed (hcp) O lattice hosting P in tetrahedral sites (PO4), Li in roughly straight lines of edgesharing octahedra (LiO6) running parallel to the b axis, and M in zigzag lines of corner-sharing octahedra (MO6) along the b axis, with Li and M in alternating b−c planes, respectively. Geometry optimizations of bulk (stoichiometric) LiMPO4 based on the empirical potential in eq 1 and DFT were performed on a unit cell of 28 atoms (4 formula units), allowing both cell shape/volume and ionic positions to relax. All optimizations were symmetry constrained. The former were carried out by means of the Newton−Raphson approach, with the Coulomb part in eq 1 computed by the Ewald summation method, and the cutoff for the remaining short-range terms set to 15 Å. For the latter, we used a 5 × 8 × 11 Gamma-centered

From ref 35.

tuning of these parameter values for LiMPO4 was attempted in our calculations, and thus this empirical potential automatically satisfies requirement c for this class of materials. Calculations based on eq 1 were performed with the GULP39,40 program as implemented in the Materials Studio software package. The DFT calculations on LiMPO4 were carried out with the VASP code41,42 within the generalized gradient approximation (GGA), using the Perdew−Burke−Ernzerhof (PBE) exchangecorrelation functional43,44 and the projector augmented wave 11205

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Monkhorst−Pack k-point mesh51 and a plane wave energy cutoff of 676 eV. Ionic/cell relaxation was performed with the conjugate gradient algorithm with a convergence tolerance on the maximum force magnitude of 10−4 eV/Å (same numerical criterion is also applied to the stress tensor components). For all DFT calculations reported in this work, electronic relaxation was performed with the RMM-DIIS algorithm41 with a convergence tolerance on the energy change of 10−6 eV. 2.2. Diffusion Barrier Calculations. We consider Li+ diffusion in LiMPO4 through the movement of vacancies and in the dilute vacancy regime. The olivine structure in Figure 2

its interior by removal of a Li+. MD simulations were performed on the resulting LiMPO4 system of 671 atoms at four temperatures, T = 700, 800, 900, and 1000 K. A leapfrog Verlet algorithm was used for the numerical integration of the equations of motion with a time step of 1 fs. The Coulomb part in eq 1 was computed with the Ewald summation method and the cutoff for the remaining short-range terms was set to 15 Å. For a given LiMPO4 system and temperature T, isothermal− isobaric ensemble (NPT) equilibration was first carried out for 1 ns. The volume was then fixed at its relaxed value and canonical ensemble (NVT) equilibration was performed for 1 ns. Finally, an NVT production run of 4 ns was carried out on this equilibrated system. For each LiMPO4 system, the trajectories from the production runs at T = 700, 800, 900, and 1000 K were used to compute the activation energy E of Li+ diffusion along the b or c axis, by first computing the values of the Li+ diffusion coefficients D(T) at these temperatures along each direction and then using these values to extract E along this direction, according to the following two-step procedure. The Li+ diffusion coefficient D can be obtained from the MD simulation using the familiar Einstein expression in the form D=

1 d lim 6N t →∞ dt

N

∑ ⟨[ri(t ) − ri(0)]2 ⟩ i=1

(2)

+

where ri is the position vector of Li i, the angular brackets represent an ensemble average (average over time origins in the MD simulation), and additional averaging over all N Li ions in the LiMPO4 system is carried out to improve the statistical precision of the calculation. As discussed above, here we focus on Li+ diffusion along only the b and c axes in LiMPO4 (paths I and II in Figure 2, respectively), and accordingly we compute D along these directions using the one-dimensional version of eq 2, Figure 2. Olivine structure of LiMPO4 (M = Mn, Fe, Co, Ni) illustrated with a 1 × 2 × 3 supercell of bulk LiFePO4. Paths I (parallel to the b axis), II (parallel to the c axis), and III for Li+ diffusion through vacancy hopping are represented by black arrows.

D=

1 d lim 2N t →∞ dt

N

∑ ⟨[qi(t ) − qi(0)]2 ⟩ i=1

(3)

where, in accordance with the adopted convention in Figure 2, q = y or z for Li+ diffusion along the b or c axis, respectively. Equation 3 implies that D along path I or II is given by half the slope of the mean square displacement (MSD) along this path in the long delay-time limit. In order to extract the energy barrier E for Li+ diffusion in LiMPO4 from a lattice site to a nearby vacancy V′Li along the b or c axis from the diffusion coefficients D(T) along this direction, a simple relationship between D and E can be derived52 by combining the harmonic transition state theory expression, p = νe−E/kT, where p is the jump or hopping frequency over the barrier, ν is the attempt frequency (typically a characteristic vibrational frequency), and k is the Boltzmann constant, with the atomistic expression for one-dimensional diffusion D = gpa2, where g is the geometrical factor and, as defined above, a is the nearest-neighbor hopping distance along the given direction, to obtain

readily reveals three possible paths30 for Li+ migration in LiMPO4 through nearest-neighbor vacancy hopping: (I) along the Li chains parallel to the b axis or in the [010] direction (parallel to the y axis), (II) between these chains and parallel to the c axis or in the [001] direction (parallel to the z axis), and (III) between these chains and in the [101] direction. Because the hopping distance a for paths I, II, and III is on average ∼3, 4.7, and 5.7 Å, respectively, path III is expected to be both the least favorable energetically31 and the most tedious to handle with the minimum energy path (MEP) DFT approach described below. Because our present interest is in assessing the empirical potential of eq 1 as a computational tool for identifying materials with specifically favorable Li+ transport capabilities, we focus here on paths I and II only. Activation energies of Li+ diffusion along paths I and II in LiMPO4 were computed from classical molecular dynamics (MD) simulations using eq 1 and from MEP DFT calculations. We describe the former calculation first. For the MD simulations, a 2 × 3 × 4 supercell of bulk LiMPO4 (e.g., ∼20.67 Å × 18.033 Å × 18.78 Å for LiFePO4) containing 672 atoms (96 formula units) was constructed, and a single VLi ′ vacancy (in Kröger−Vink notation) was created in

D = D0 e−E / kT

(4)

where D0 = gνa2. Equation 4 can now be recast as log D = log D0 − 11206

⎛ ⎞ 103 E ⎜ ⎟ ⎝ 103k ln 10 ⎠ T

(5)

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Table 2. Lattice Parameters a, b, and c of LiMPO4 (M = Mn, Fe, Co, Ni)a eq 1 a (Å) b (Å) c (Å)

10.6018 (−1.64%) 6.1952 (−1.65%) 4.8444 (−2.27%) 10.4057 (−0.66%) 6.0830 (−1.19%) 4.8094 (−2.44%) 10.3426 (−1.40%) 6.0319 (−1.89%) 4.8105 (−2.57%) 10.2459 (−2.18%) 5.9876 (−2.29%) 4.7942 (−2.52%)

DFT LiMnPO4 10.5537 (−1.18%) 6.1476 (−0.87%) 4.7800 (−0.92%) LiFePO4 10.4011 (−0.61%) 6.0434 (−0.53%) 4.7399 (−0.96%) LiCoPO4 10.3174 (−1.15%) 5.9279 (−0.13%) 4.7393 (−1.05%) LiNiPO4 10.1286 (−1.01%) 5.9211 (−1.15%) 4.7204 (−0.94%)

experiment

other

10.431(9)b 6.0947(6)b 4.7366(5)b

10.5401c 6.0874c 4.6878c

10.3377(5)d 6.0112(2)d 4.6950(2)d

10.3713,c,e 10.362f 6.0216,c,e 5.983f 4.6695,c,e 4.680f

10.2001(6)g 5.9199(4)g 4.690(2)g

10.2428c 5.9093c 4.6418c

10.0275(3)b 5.8537(2)b 4.6763(2)b

10.1353c 5.8432c 4.6257c

a

Values obtained from our geometry optimizations based on the empirical potential of eq 1 and DFT, from experiment, and from other computational studies employing empirical potentials customized for LiMPO4. Relative errors with respect to experimental values are shown in parentheses. bFrom ref 50. cFrom ref 31. dFrom ref 47. eFrom ref 23. fFrom ref 24. gFrom ref 49.

the VTST interface54 to VASP, to determine the MEP and associated activation energy barrier E for Li+ diffusion in LiMPO4 from a lattice site to a nearby vacancy V′Li (created by removal of a Li+ from its lattice site) along the b or c axis (paths I or II in Figure 2, respectively). To compute E for Li+ diffusion along a given axis, starting with two identical 1 × 2 × 3 supercells of bulk LiMPO4 optimized as described above, a V′Li vacancy is introduced at an appropriate location in each, and geometry optimization is performed on the resulting supercells of 167 atoms, with the supercell parameters fixed to their optimized bulk values but allowing all ions to relax, to give the two end-point configurations (end-point images) required at the start of a CI-NEB calculation. This geometry optimization was carried out at the Gamma point with a plane wave energy cutoff of 400 eV. Ionic relaxation was performed with the conjugate gradient algorithm with a convergence tolerance on the maximum force magnitude of 10−4 eV/Å. Eight initial intermediate configurations (images) were generated by linear interpolation between the end-point images. Relaxation of the intermediate images to the MEP was performed at the Gamma point with a plane wave energy cutoff of 400 eV (same as above for the endpoint images), using the quick-min algorithm55 with a convergence tolerance on the maximum force magnitude of 10−4 eV/Å. Asymptotic convergence of E was almost always reached before this rather stringent condition was satisfied, at which point the relaxation was stopped.

implying that, for D along a desired direction, E along this direction can be extracted from the slope of the Arrhenius plot of log D vs 103/T. The geometrical factor g is a constant that depends on the lattice type and accounts for (i.e., is the reciprocal of) the number of equivalent nearest-neighbor jumps of the vacancy. It is typically of order unity and in the present case, where vacancy jumps are constrained to only one dimension along either path I or II, g = 1/2. In summary then, for each LiMPO4 system, we first compute the MSD of Li+ along the b or c axis (averaged over all Li ions in the system) as a function of the delay-time over the entire production phase (4000 ps) of the simulations at T = 700, 800, 900, and 1000 K, with successive time origins separated by a single sampling time step of size 0.4 ps. We have experimented with various sampling time step sizes to ensure that the sampling resolution is sufficiently accurate. We neglect the trailing part of the MSD vs delay-time curve (which typically suffers from poor statistical precision because of the progressively decreasing sample size) spanning the last 3700 ps and perform an unweighted linear regression fit on the first 300 ps of the curve, which gives D(T) along the desired axis as half the slope of the fitted line, as per eq 3. We then perform an unweighted linear regression fit on the Arrhenius plot of the four (log D, 103/T) points (corresponding to T = 700, 800, 900, 1000 K) and extract the energy barrier E along the desired axis from the slope of the fitted line according to eq 5. The LiMPO4 intermediate results from the above procedure are all reported in section 3. For the MEP DFT calculations of the activation energies of Li+ diffusion along paths I and II in LiMPO4, a 1 × 2 × 3 supercell of bulk LiMPO4 (e.g., ∼10.33 Å × 12.022 Å × 14.085 Å for LiFePO4) as in Figure 2, containing 168 atoms (24 formula units) was constructed, and geometry optimization was performed on it, using a 4 × 3 × 3 Gamma-centered Monkhorst−Pack k-point mesh51 and a plane wave energy cutoff of 520 eV and allowing both cell parameters and ionic positions to relax. Ionic/cell relaxation was performed in the same fashion as described for the LiMPO4 unit cell in section 2.1. Following optimization, we used the climbing image nudged elastic band (CI-NEB) method,53 as implemented in

3. RESULTS AND DISCUSSION We first compare the ground state structures of bulk LiMPO4 (M = Mn, Fe, Co, Ni) obtained from the geometry optimizations with the potential of eq 1 to those generated by the geometry optimizations with DFT and to experimental structures. We then validate the Li+ diffusion barriers in LiMPO4 obtained from our DFT-based CI-NEB calculations against those from other DFT studies and some experimental results and subsequently compare the barriers obtained from the classical MD simulations based on eq 1 with these validated DFT values. 11207

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Figure 3. Activation energy E of Li+ diffusion along the b axis in LiMPO4 (M = Mn, Fe, Co, Ni) obtained from our DFT calculations with the CI-NEB method. The two end-point minima images and eight intermediate images used in the calculations are represented by dots.

Figure 4. Same as Figure 3 but for Li+ diffusion along the c axis in LiMPO4 (M = Mn, Fe, Co, Ni).

empirical potential of eq 1. These values are in good overall agreement with the DFT and experimental ones, with their relative errors with respect to experiment (parenthesized quantities in Table 2) generally close to but, as expected from the empirical nature of the potential and the strict satisfaction of requirement c for LiMPO4, larger in magnitude than the relative errors of the DFT values with respect to experiment. The good agreement between the geometries optimized with the empirical potential of eq 1 and those from experiment and DFT reflects the ability of this potential to reproduce accurately the forces in the neighborhood of the potential energy surface minimum.

3.1. Ground State Structures. Table 2 contains the values of the lattice parameters of LiMPO4 obtained from our DFT based geometry optimizations, performed as described in section 2.1. Experimental values of the lattice parameters are also reported in Table 2. We note that our DFT values of the lattice parameters are in good agreement with the experimental values and only slightly larger throughout, as expected from our usage of a GGA exchange-correlation functional (PBE43,44). Table 2 contains also the values of the lattice parameters obtained from our geometry optimizations based on the 11208

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Table 3. Computed Activation Energies of Li+ Diffusion along the b and c Axes (Paths I and II in Figure 2, respectively) in LiMPO4 (M = Mn, Fe, Co, Ni)a M

E* (eV)

E† (eV)

Mn Fe Co Ni

0.264 0.228 0.539 0.473

0.315 0.321 0.283 0.225

Mn Fe Co Ni

2.029 2.616 2.384 2.025

2.072 2.301 2.425 2.578

E** (eV)

Along the b axis 0.62b 0.55,b,d 0.57,e 0.420f 0.49b 0.44b Along the c axis 2.83b 2.89b,d 3.28b 3.49b

E†† (eV) 0.250c 0.270,c 0.32g 0.360c 0.130c

2.5,c 2.27g

E* from our MD simulations using the empirical potential of eq 1, E† from our DFT-based CI-NEB calculations, E** from other studies with custom parametrized empirical potentials, and E†† from other DFT studies where available. bFrom ref 31. cFrom ref 30. dFrom ref 23. eFrom ref 25. fFrom ref 24. gFrom ref 27. a

eq 1 (which has not been customized for LiMPO4, as emphasized early in section 2) and the experimental lattice parameters is by comparison quite good. With this partial validation of the empirical potential of eq 1 and its parameters in Table 1, based on the above LiMPO4 structural results, we now proceed to the more crucial question concerning the ability of this potential to predict at least qualitatively correct Li+ diffusion barriers. 3.2. Activation Energies. Figures 3 and 4 show the energy profiles for Li+ diffusion along the b and c axes in LiMPO4 (paths I and II in Figure 2), respectively, from our DFT-based CI-NEB calculations detailed in section 2.2. The figures also display the relative energies of the eight intermediate images along the MEP, including the numerical value of the energy barrier E. These results imply that, in the absence of extrinsic or other intrinsic defects, Li+ diffusion through vacancy hopping in LiMPO4 occurs predominantly along the one-dimensional Li channels parallel to the b axis (i.e., along path I), with a relatively high energy barrier inhibiting Li+ migration between these channels (i.e., along path II). Figure 5 shows, for each LiMPO4 system, a view from the top along the c axis of an overlay of the ten images featured in the corresponding energy profile from Figure 3. As illustrated in this figure, we have observed that the actual Li+ migration trajectory between adjacent sites along the b axis is bent rather than straight, for all four LiMPO4 systems, resulting in a zigzag Li+ migration along path I. A close-up look reveals also the slight distortion of the local atomic environment of the migrating Li+, as evidenced by small relative displacements of corresponding atoms in this environment from the different images that would otherwise overlap perfectly (not all such displacements are discernible due to the reduced scale of the figure). All of these conclusions are in general agreement with results from previous computational studies23,27,30,31 and experiment,56−59 although some experimental results suggesting mainly two-dimensional Li+ conductivity in the b−c plane have also been reported.60,61 Table 3 compares the activation energies of Li+ diffusion along the b and c axes in LiMPO4, produced by our DFT-based CI-NEB calculations (E†) with those from other DFT calculations (E††). We note that the activation energies E† of Li+ diffusion along the c axis in LiMnPO4, LiCoPO4, and LiNiPO4 represent new DFT results, which serve as a reference

Figure 5. Top view along the c axis (z axis) of a superposition of the ten configurations or images (two images representing the end-point minima and eight intermediate images) along the MEP from our DFT CI-NEB calculations of Li+ diffusion barriers along the b axis in LiMPO4 (M = Mn, Fe, Co, Ni). The Li atoms are depicted larger than the M, P, and O atoms, for clarity. The migration path of the Li+ is clearly seen to be bent, rather than straight, in all cases. Note that the path is not in the plane of the figure.

Table 2 also includes lattice parameters of LiMPO 4 computed in other investigations using system-tuned empirical potentials. Given that the parameters of such empirical potentials are by definition designed to produce high quality agreement between the predicted and experimental lattice parameters (among other properties) of specifically LiMPO4, as verified from Table 2, the agreement between the lattice parameters predicted on the basis of the empirical potential in 11209

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Figure 6. Final atomic configurations from the canonical ensemble MD simulations carried out on LiMPO4 (M = Mn, Fe, Co, Ni) with the empirical potential of eq 1 at temperature T = 800 K. The box represents the 2 × 3× 4 supercell of LiMPO4 making up the MD system (see section 2.2). The Li atoms are depicted as dots, and the M, P, and O atoms (confined to the box) are suppressed for clarity. The Li+ diffusion is observed to be predominantly parallel to the b axis in all four systems.

for judging the accuracy of not just the corresponding values we obtained using the empirical potential of Pedone et al., as carried out below, but also those obtained by others using different empirical potentials. Taking into account differences in computational parameter values (e.g., cell size, k-point mesh, etc., see section 2) between our work and the other calculations, the agreement between E† and E†† is overall good. In particular, there is very good agreement with the careful calculations of Hoang and Johannes27 for LiFePO4 (0.321 vs 0.32 along path I, and 2.301 vs 2.27 along path II), which were carried out with antiferromagnetic spin polarization and using the GGA+U approach. Therefore, as anticipated early in section 2, our GGA/DFT activation energies obtained with ferromagnetic spin polarization are sufficiently accurate to serve as a reference with which we compare below the activation energies generated with the empirical potential of eq 1. Figure 6 provides snapshots of the Li atomic configurations taken at the end of the 4 ns production runs of the canonical ensemble MD simulations at T = 800 K, performed on LiMPO4

using the empirical potential of eq 1, as detailed in section 2.2. From this figure, the Li+ diffusion in LiMPO4 is clearly seen to occur mainly along the b axis. On a more quantitative level, Figures 7−10 show the MSD of Li+ along the b axis in LiMPO4 computed from the MD simulations at temperatures T = 700, 800, 900, and 1000 K, together with their linear regression fits. The insets show Arrhenius plots for the Li+ diffusion coefficients obtained from the slopes of these fitted lines, together with their linear regression fits and energy barriers given by E = −103 k ln 10 × (slope of fitted line) according to eq 5. Similar remarks apply to Figures 11−14, but for Li+ diffusion along the c axis in LiMPO4. The energy barrier values obtained from these MD simulations again favor Li+ diffusion along the b axis with relatively negligible cross channel migration. Thus, the picture of Li+ transport emerging from our MD simulations based on eq 1, as occurring primarily along the b axis in LiMPO4, is in agreement with the findings from the above DFT-based CI-NEB results. 11210

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Figure 7. Activation energy E of Li+ diffusion along the b axis in LiMnPO4 calculated from our canonical ensemble MD simulations based on the empirical potential of eq 1 at temperatures T = 700, 800, 900, and 1000 K. Solid lines represent the MSD of Li+ along the b axis computed from the simulations and dashed lines their unweighted linear regression fits. The values of the Li+ diffusion coefficients D(T) along the b axis obtained from the slopes of the latter are used to generate the dots in the Arrhenius plot inset. The dashed line there represents their unweighted linear regression fit and E is extracted from its slope according to eq 5.

Figure 10. Same as Figure 7 but for LiNiPO4.

Figure 11. Same as Figure 7 but for Li+ diffusion along the c axis in LiMnPO4. The MSD line for T = 700 K showed no discernible slope and accordingly was not used in generating the inset fit.

Figure 8. Same as Figure 7 but for LiFePO4.

Figure 12. Same as Figure 11 but for LiFePO4.

Furthermore, Table 3 provides a comparison of the numerical values of the activation energies of Li+ diffusion in LiMPO4 in the dilute vacancy limit along the b and c axes, computed from the MD simulations (E*) and the DFT CINEB approach (E†). The values of E* and E† are seen to be in good overall agreement and for some cases in very good

Figure 9. Same as Figure 7 but for LiCoPO4. 11211

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

they did, for example, for the lattice parameter values in Table 2. The issue of agreement in trend between E* and E† is more subtle. Whenever the DFT barrier values along some direction are close to each other, as is the case here for LiMPO4, two complications arise. First, consider from Table 3 the trends of the activation energies of Li+ diffusion along the b axis in LiMPO4 for E† (Ni < Co < Mn < Fe), E†† (Ni < Mn < Fe < Co), and E** (Ni < Co < Fe < Mn) from ref 31. Which of these three different trends should be taken as reference to compare the trend in E* against? The small differences in barrier values have therefore introduced an uncertainty in the accurately computed trends. Second, even if a unique reference trend in barrier values could be identified, the required qualitative agreement in values between E* and E† is unlikely to also yield an agreement in trend between them, again because of the proximity of E† values to each other. Fortunately, these complications are of no consequence for our desired computational tool. For a tool intended for scanning large numbers of cathode materials and selecting those with desirable Li+ transport capabilities (e.g., low activation barriers) for more sophisticated DFT calculations, the relevant feature from Table 3 is that the empirical potential of eq 1 predicts correctly the observed order of magnitude difference between the DFT activation energies along the b and c axes. The ability to predict the correct ordering (sometimes unknown as shown above) of closely spaced activation barriers along some direction serves no purpose in such a computational tool. The quality of agreement between E* and E† in Table 3 is sufficiently good to characterize the values of E* as at least qualitatively correct compared with E† and satisfy requirement a for the present LiMPO4 test cases, as set out in section 1, leading to a positive initial assessment of the empirical potential of eq 1 as a rapid portable computational tool for singling out materials with favorable Li+ transport capabilities. With this favorable initial assessment in hand, we are currently expanding the range of test cathode materials for this potential toward the more comprehensive assessment mentioned in section 1. Finally, Table 4 offers a comparison between the diffusion coefficients of Li+ along the b axis in LiMPO4 at T = 1000 K, computed with the DFT activation energies E† along the b axis from Table 3 by means of eq 4 and those obtained directly from the MD simulations based on eq 1 at T = 1000 K (see Figures 7−10) as detailed in section 2.2. For this purpose, we took ν = 1012 Hz and a = 3 Å (average Li−Li distance along the Li chain parallel to the b axis) in eq 4. Again, Table 4 shows overall good agreement between the two sets of values, at least good enough to satisfy requirement a for LiMPO4, our core objective here.

Figure 13. Same as Figure 11 but for LiCoPO4.

Figure 14. Same as Figure 11 but for LiNiPO4.

Table 4. Diffusion Coefficients D of Li+ along the b Axis (Path I in Figure 2) in LiMPO4 (M = Mn, Fe, Co, Ni) at Temperature T = 1000 Ka M

DMD (×10−5 cm2/s)

DDFT (×10−5 cm2/s)

Mn Fe Co Ni

0.91 1.53 1.28 1.12

1.16 1.08 1.68 3.30

a

DMD is obtained directly from our MD simulations based on eq 1 at T = 1000 K (see Figures 7−10) and DDFT is computed using E† along the b axis from Table 3 in eq 4.

4. CONCLUSION We have reported ground state structure and Li+ diffusion barrier calculations on the phosphate olivines, LiMPO4 (M = Mn, Fe, Co, Ni), which are currently attracting a great deal of research interest as potential next generation Li-ion battery cathode materials. The calculations were carried out both classically and ab initio. The classical calculations were based on the empirical potential of Pedone et al.35 and, for computing the activation energies for Li+ diffusion, MD simulations were performed in the canonical ensemble at four different temperatures for each LiMPO4 system, allowing for the computation of the Li+ diffusion coefficients at these temper-

agreement. Achieving exact agreement between E* and E† is simply out of the question for any empirical potential and was never our aim to begin with. Nevertheless, we observe from Table 3 that the quality of agreement between E† and E* is in most cases higher than that between E† and E**, where E** is obtained from other classical studies using various empirical potentials custom fitted for LiMPO4. This observation is somewhat surprising given that E* was generated here by the empirical potential of eq 1, which strictly satisfied requirement c for LiMPO4 with the original parameter values in Table 1 unaltered throughout, whereas E** was produced by empirical potentials tailored to give more accurate results for LiMPO4, as 11212

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

(8) Wen, C. J.; Huggins, R. A. Thermodynamic Study of the LithiumTin System. J. Electrochem. Soc. 1981, 128, 1181−1187. (9) Wen, C. J.; Huggins, R. A. Electrochemical Investigation of the Lithium-Gallium System. J. Electrochem. Soc. 1981, 128, 1636−1641. (10) Rao, B. M. L.; Francis, R. W.; Christopher, H. A. LithiumAluminum Electrode. J. Electrochem. Soc. 1977, 124, 1490−1492. (11) Murphy, D. W.; Disalvo, F. J.; Carides, J. N.; Waszczak, J. V. Topochemical Reactions of Rutile Related Structures with Lithium. Mater. Res. Bull. 1978, 13, 1395−1402. (12) Murphy, D. W.; Carides, J. N. Low Voltage Behavior of Lithium/Metal Dichalcogenide Topochemical Cells. J. Electrochem. Soc. 1979, 126, 349−351. (13) Lazzari, M.; Scrosati, B. Cyclable Lithium Organic Electrolyte Cell Based on Two Intercalation Electrodes. J. Electrochem. Soc. 1980, 127, 773−774. (14) Tarascon, J. M.; Armand, M. Issues and Challenges Facing Rechargeable Lithium Batteries. Nature 2001, 414, 359−367. (15) Mizushima, K.; Jones, P. C.; Wiseman, P. J.; Goodenough, J. B. LixCoO2 (0 < x ≤ 1): A New Cathode Material for Batteries of High Energy Density. Solid State Ionics 1981, 3−4, 171−174. (16) Mizushima, K.; Jones, P. C.; Wiseman, P. J.; Goodenough, J. B. LixCoO2 (0 < x ≤ 1): A New Cathode Material for Batteries of High Energy Density. Mater. Res. Bull. 1980, 15, 783−789. (17) Yazami, R.; Touzain, P. A Reversible Graphite Lithium Negative Electrode for Electrochemical Generators. J. Power Sources 1983, 9, 365−371. (18) Scrosati, B. History of Lithium Batteries. J. Solid State Electrochem. 2011, 15, 1623−1630. (19) Hayner, C. M.; Zhao, X.; Kung, H. H. Materials for Rechargeable Lithium-Ion Batteries. Annu. Rev. Chem. Biomol. Eng. 2012, 3, 445−471. (20) Shi, Z.; Liu, H.; Zhang, J. In Lithium Batteries: Research, Technology and Applications; Dahlin, G. R., Strøm, K. E., Eds.; Electrical Engineering Developments Series; Nova Science Publishers: New York, 2010. (21) Padhi, A. K.; Nanjundaswamy, K. S.; Goodenough, J. B. Phospho-olivines as Positive-Electrode Materials for Rechargeable Lithium Batteries. J. Electrochem. Soc. 1997, 144, 1188−1194. (22) Chung, S. Y.; Bloking, J. T.; Chiang, Y. M. Electronically Conductive Phospho-olivines as Lithium Storage Electrodes. Nat. Mater. 2002, 1, 123−128. (23) Islam, M.; Driscoll, D.; Fisher, C.; Slater, P. Atomic-Scale Investigation of Defects, Dopants, and Lithium Transport in the LiFePO4 Olivine-Type Battery Material. Chem. Mater. 2005, 17, 5085−5092. (24) Kuss, C.; Liang, G.; Schougaard, S. B. Atomistic Modeling of Site Exchange Defects in Lithium Iron Phosphate and Iron Phosphate. J. Mater. Chem. 2012, 22, 24889−24893. (25) Adams, S.; Rao, R. P. Simulated Defect and Interface Engineering for High Power Li Electrode Materials. Solid State Ionics 2011, 184, 57−61. (26) Boulfelfel, S. E.; Seifert, G.; Leoni, S. Atomistic Investigation of Li+ Diffusion Pathways in the Olivine LiFePO4 Cathode Material. J. Mater. Chem. 2011, 21, 16365−16372. (27) Hoang, K.; Johannes, M. Tailoring Native Defects in LiFePO4: Insights from First-Principles Calculations. Chem. Mater. 2011, 23, 3003−3013. (28) Ouyang, C. Y.; Shi, S. Q.; Wang, Z. X.; Huang, X. J.; Chen, L. Q. First-Principles Study of Li Ion Diffusion in LiFePO4. Phys. Rev. B 2004, 69, No. 104303. (29) Yang, J.; Tse, J. S. Li Ion Diffusion Mechanisms in LiFePO4: An Ab Initio Molecular Dynamics Study. J. Phys. Chem. A 2011, 115, 13045−13049. (30) Morgan, D.; Van der Ven, A.; Ceder, G. Li Conductivity in LixMPO4 (M = Mn, Fe, Co, Ni) Olivine Materials. Electrochem. Solid State Lett. 2004, 7, A30−A32. (31) Fisher, C. A. J.; Prieto, V. M. H.; Islam, M. S. Lithium Battery Materials LiMPO4 (M = Mn, Fe, Co, and Ni): Insights into Defect

atures along the b and c axes and extraction of the desired activation barrier along each axis. The first-principles calculations were based on DFT, and the CI-NEB method53 was employed for determining the MEP and corresponding energy barrier for Li+ diffusion through vacancy hopping along the b and c axes in each LiMPO4 system. We showed that the ground state structures of LiMPO4 from our geometry optimizations with the potential of eq 1 are in good agreement with the experimental structures and our DFT ones. After validating our DFT barrier results by comparison with experiment and previous DFT investigations of Li+ transport in LiMPO4, we showed that our classical energy barriers are in fairly good agreement with our DFT values for the LiMPO4 test systems considered here. Our broad objective was to provide an initial assessment, based on the LiMPO4 test cases, of whether the empirical potential of Pedone et al. could potentially provide a fast, portable, and at least qualitatively accurate computational tool for picking out, from many candidates, potential cathode materials with favorable Li+ transport capabilities, which subsequently could be examined with greater rigor but also at greater computational expense with DFT studies similar to the ones we performed and reported here. Such a tool can be exploited to reduce wasteful expensive computations, by focusing CPU-intensive DFT studies on those materials with the more promising Li+ transport capabilities. Given that for LiMPO4 the agreement between our classical and DFT Li+ diffusion energy barriers is at least qualitatively good, indeed typically better, as we have shown, than that between the latter and barrier values generated by other researchers with custom parametrized empirical potentials, our initial assessment of the empirical potential of Pedone et al. is promising. Building on these early encouraging LiMPO4 results, we are currently performing and will report similar calculations on additional cathode materials with the aim of providing a more broadly based assessment of this empirical potential.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] or [email protected]. Phone: +61 (0)2 9717 9963. Notes

The authors declare no competing financial interest.



REFERENCES

(1) Armand, M.; Tarascon, J. M. Building Better Batteries. Nature 2008, 451, 652−657. (2) Nagaura, T.; Tozawa, K. Lithium Ion Rechargeable Battery. Prog. Batteries Solar Cells 1990, 9, 209−217. (3) Ozawa, K. Lithium-Ion Rechargeable Batteries with LiCoO2 and Carbon Electrodes - The LiCoO2/C System. Solid State Ionics 1994, 69, 212−221. (4) Whittingham, M. S. Electrical Energy Storage and Intercalation Chemistry. Science 1976, 192, 1126−1127. (5) Weppner, W.; Huggins, R. A. Determination of the Kinetic Parameters of Mixed-Conducting Electrodes and Application to the System Li3Sb. J. Electrochem. Soc. 1977, 124, 1569−1578. (6) Wen, C. J.; Boukamp, B. A.; Huggins, R. A.; Weppner, W. Thermodynamic and Mass Transport Properties of “LiAl”. J. Electrochem. Soc. 1979, 126, 2258−2266. (7) Wen, C. J.; Huggins, R. A. Chemical Diffusion in Intermediate Phases in the Lithium-Silicon System. J. Solid State Chem. 1981, 37, 271−278. 11213

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214

The Journal of Physical Chemistry C

Article

Association, Transport Mechanisms, and Doping Behavior. Chem. Mater. 2008, 20, 5907−5915. (32) Shoichet, B. K. Virtual Screening of Chemical Libraries. Nature 2004, 432, 862−865. (33) Jain, A.; Hautier, G.; Moore, C. J.; Ong, S. P.; Fischer, C. C.; Mueller, T.; Persson, K. A.; Ceder, G. A High-Throughput Infrastructure for Density Functional Theory Calculations. Comput. Mater. Sci. 2011, 50, 2295−2310. (34) Curtarolo, S.; Hart, G. L. W.; Nardelli, M. B.; Mingo, N.; Sanvito, S.; Levy, O. The High-Throughput Highway to Computational Materials Design. Nat. Mater. 2013, 12, 191−201. (35) Pedone, A.; Malavasi, G.; Menziani, M. C.; Cormack, A. N.; Segre, U. A New Self-Consistent Empirical Interatomic Potential Model for Oxides, Silicates, and Silica-Based Glasses. J. Phys. Chem. B 2006, 110, 11780−11795. (36) Afify, N. D.; Mountjoy, G. Molecular-Dynamics Modeling of Eu3+-Ion Clustering in SiO2 Glass. Phys. Rev. B 2009, 79, No. 024202. (37) Pedone, A. Properties Calculations of Silica-Based Glasses by Atomistic Simulations Techniques: A Review. J. Phys. Chem. C 2009, 113, 20773−20784. (38) Ori, G.; Montorsi, M.; Pedone, A.; Siligardi, C. Insight into the Structure of Vanadium Containing Glasses: A Molecular Dynamics Study. J. Non-Cryst. Solids 2011, 357, 2571−2579. (39) Gale, J. D. GULP: A Computer Program for the SymmetryAdapted Simulation of Solids. J. Chem. Soc., Faraday Trans. 1997, 93, 629−637. (40) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simul. 2003, 29, 291−341. (41) Kresse, G.; Furthmuller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (42) Kresse, G.; Furthmuller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169. (43) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple (vol 77, pg 3865, 1996). Phys. Rev. Lett. 1997, 78, 1396. (45) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953. (46) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758. (47) Rousse, G.; Rodriguez-Carvajal, J.; Patoux, S.; Masquelier, C. Magnetic Structures of the Triphylite LiFePO4 and of its Delithiated Form FePO4. Chem. Mater. 2003, 15, 4082−4090. (48) Ong, S. P.; Chevrier, V. L.; Hautier, G.; Jain, A.; Moore, C.; Kim, S.; Ma, X.; Ceder, G. Voltage, Stability and Diffusion Barrier Differences Between Sodium-Ion and Lithium-Ion Intercalation Materials. Energy Environ. Sci. 2011, 4, 3680−3688. (49) Kubel, F. Crystal-Structure of Lithium Cobalt Double Orthophosphate, LiCoPO4. Z. Kristallogr. 1994, 209, 755. (50) Garcia-Moreno, O.; Alvarez-Vega, M.; Garcia-Alvarado, F.; Garcia-Jaca, J.; Gallardo-Amores, J. M.; Sanjuan, M. L.; Amador, U. Influence of the Structure on the Electrochemical Performance of Lithium Transition Metal Phosphates as Cathodic Materials in Rechargeable Lithium Batteries: A New High-Pressure Form of LiMPO4 (M = Fe and Ni). Chem. Mater. 2001, 13, 1570−1576. (51) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. (52) Kittel, C. Introduction to Solid State Physics, 6th ed.; Wiley: New York, 1986. (53) Henkelman, G.; Uberuaga, B. P.; Jonsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (54) See http://theory.cm.utexas.edu/vtsttools/. (55) Sheppard, D.; Terrell, R.; Henkelman, G. Optimization Methods for Finding Minimum Energy Paths. J. Chem. Phys. 2008, 128, No. 134106.

(56) Nishimura, S.-i.; Kobayashi, G.; Ohoyama, K.; Kanno, R.; Yashima, M.; Yamada, A. Experimental Visualization of Lithium Diffusion in LixFePO4. Nat. Mater. 2008, 7, 707−711. (57) Li, J.; Yao, W.; Martin, S.; Vaknin, D. Lithium Ion Conductivity in Single Crystal LiFePO4. Solid State Ionics 2008, 179, 2016−2019. (58) Baker, P. J.; Franke, I.; Pratt, F. L.; Lancaster, T.; Prabhakaran, D.; Hayes, W.; Blundell, S. J. Probing Magnetic Order in LiMPO4 (M = Ni, Co, Fe) and Lithium Diffusion in LixFePO4. Phys. Rev. B 2011, 84, No. 174403. (59) Sugiyama, J.; Nozaki, H.; Harada, M.; Kamazawa, K.; Ikedo, Y.; Miyake, Y.; Ofer, O.; Mansson, M.; Ansaldo, E. J.; Chow, K. H.; et al. Diffusive Behavior in LiMPO4 with M = Fe, Co, Ni Probed by MuonSpin Relaxation. Phys. Rev. B 2012, 85, No. 054111. (60) Amin, R.; Maier, J. Effect of Annealing on Transport Properties of LiFePO4: Towards a Defect Chemical Model. Solid State Ionics 2008, 178, 1831−1836. (61) Amin, R.; Maier, J.; Balaya, P.; Chen, D. P.; Lin, C. T. Ionic and Electronic Transport in Single Crystalline LiFePO4 Grown by Optical Floating Zone Technique. Solid State Ionics 2008, 179, 1683−1687.

11214

dx.doi.org/10.1021/jp5004402 | J. Phys. Chem. C 2014, 118, 11203−11214