Cooperative Mechanism for the Diffusion of Li+ Ions in LiMgSO4F

Jul 18, 2012 - School of Chemistry, Trinity College Dublin, Dublin 2, Ireland .... diffusion of Li+ ions, even at high temperatures.46,47 For this rea...
1 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCC

Cooperative Mechanism for the Diffusion of Li+ Ions in LiMgSO4F Mathieu Salanne,*,† Dario Marrocchelli,*,‡ and Graeme W. Watson‡ †

UPMC Univ-Paris06 and CNRS, UMR 7195, PECSA, F-75005 Paris, France, and Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, France ‡ School of Chemistry, Trinity College Dublin, Dublin 2, Ireland S Supporting Information *

ABSTRACT: Extensive molecular dynamics simulations are performed on tavorite-structured LiMgSO4F, in order to analyze the diffusion mechanism of the Li+ ions. In a first step, an interaction potential, including polarization effects, is parametrized from density functional theory calculations. This is then tested by reproducing experimental properties of the material, such as structural parameters (lattice constants and interatomic distances) and ionic conductivity to a high degree of accuracy. Next, the conduction mechanism is studied: the diffusion of the Li+ ions occurs via correlated hops inside diffusion channels, which reflects in the much larger values calculated for their collective diffusion coefficient compared to the individual ones (by 1 order of magnitude on average). The consequences (both scientific and technical) of this finding are discussed.

1. INTRODUCTION Spurred by the enormous needs in electricity storage devices, the search of new materials for lithium-ion batteries has received much attention in the recent years.1 Among the noticeable advances, the introduction of cathode materials based on polyoxoanionic species such as phosphates (e.g., LiFePO4),2−4 pyrophosphates,5 or, more recently, sulphates (LiFeSO4F)6,7 seems to be one of the most promising. The interest in battery electrodes has gone much beyond the classical electrochemistry field, and numerous developments have been made for characterizing the materials. Nowadays, the vast array of experimental methods encompasses many diffraction or spectroscopy techniques such as X-ray scattering,8,9 Mossbauer,10 or NMR.11−13 Computational methods have also made a very useful contribution to the field. In particular, systematic first-principles calculations performed in the framework of the density functional theory (DFT) have increasingly been used to identify cathode materials with optimal properties.14,15 The always rising computational power available and the design of powerful algorithms have allowed systematic studies on thousands of compounds; recent illustrations of the capabilities of this approach are the evaluations of phosphate-based compounds16 and tavorite-structured materials17 using high-throughput calculations. DFT calculations have also been used in order to interpret NMR spectroscopy experiments.18,19 An elegant approach, complementary to DFT calculations, is the use of a qualitative picture of electronic band structure, i.e., molecular orbital diagrams, in order to correlate the nature of the electronic states involved in redox reactions with the capacity of the materials.20 Such an analysis was also used to rationalize the metastability of the tavorite phase of LiMnSO4F.21 © 2012 American Chemical Society

Although DFT calculations have shown an unprecedented predictive ability concerning static properties, such as the cathode materials voltage, the situation is somewhat more ambiguous for dynamic properties such as the diffusion coefficient of Li+ ions and the electrical conductivity. In fact, diffusion coefficients D are usually estimated from D = a 2 ν* exp( −Eact /kBT )

(1)

where a is the Li+ ion hopping length, ν* is the hop attempt frequency, and Eact is the activation barrier. Such an estimation relies on several hypothesis, among which the most stringent is the assumption that Li+ ions are not interacting with each other.22 It is generally inferred that this hypothesis is valid for materials that are either fully delithiated or entirely lithiated, i.e., in the conditions of absence of vacancies.23 Among the three parameters of eq 1, one only (Eact) is effectively calculated by DFT. Indeed, a is chosen based on the geometry of the materials, and ν* is usually taken as equal to 1012 Hz. Unfortunately, this approach seems to provide diffusion coefficients that are too high compared to experiments. In LiFePO4, for example, it provides a value of around 10−8 cm2 s−1 at room temperature23 (or around 10−6 cm2 s−1 at 200 °C), which is a couple of orders of magnitude higher than the reported experimental value (which is of 10−8 cm2 s−1 at 200 °C24 for single crystals). It is worth noting that the predicted diffusion coefficient for Li+ ions is only 1 order of magnitude smaller than the typical diffusion coefficients in ionic liquids at room temperature25 and that such values imply that the Received: May 16, 2012 Revised: July 12, 2012 Published: July 18, 2012 18618

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C

Article

from the computational point of view, and the accessible simulation times for our systems of interest does not go beyond a few tenths of a picosecond, which is not sufficient to observe diffusion of Li+ ions, even at high temperatures.46,47 For this reason, in recent years, some physically motivated models for ionic materials have been introduced.48 Several levels of complexity can be used for the analytical form, depending on the importance of multipolar effects or of ion distorsions. In the present work, the interaction potential is decomposed into four terms:

diffusive régime is reached by the particles on the time scales of 10−100 ns. The difference between the calculated and experimental values was explained based on a model accounting for the presence of blocking sites along the one-dimensional channel paths.26 Classical molecular dynamics (MD) simulations should, in principle, provide an alternative route to study the dynamics in cathode materials. This method yields the trajectory of the atoms from which the diffusion coefficient can be calculated from the mean-squared displacements D = lim

t →∞

1 ⟨|δ ri(t )|2 ⟩ 6t

Vtotal = Vcharge + Vdispersion + Vrepulsion + Vinduction (2)

(3)

The first term is the direct Coulomb interaction between two ions i and j

provided that the diffusive régime has been reached. In this expression, δri(t) is the displacement of a given ion at time t and the brackets denote an ensemble average. Here, the calculation of D is direct in the sense that it does not rely on any model. Another advantage of classical MD simulations is that they allow the study of relatively big systems (thousands of atoms) for long time scales (several nanoseconds). This approach has successfully been used in studies of numerous superionic materials such as zirconia- or ceria-based electrolytes27−29 and novel solid oxide fuel cell cathode materials.30−32 Several cathode materials for Liion batteries have therefore recently been characterized by MD simulations.33−36 In this work, we perform DFT calculations on the tavorite structure of LiMgSO4F and use them to parametrize a polarizable interaction potential. We then perform extensive MD simulations of this material, involving large simulation cells. LiMgSO4F was selected because it is a structural analogue of LiFeSO4F, a promising material for future battery applications, that, however, does not exhibit electronic conductivity. This allows us to focus on the ionic conduction properties only. We expect the conduction mechanism of this material to be fairly similar to that of stoichiometric LiFeSO4F, i.e., when Fe has the same valency as Mg2+. This methodology, in principle, does not allow the study of Li vacancies in this material, though possible developments in this direction are discussed at the end of the article. We test our potential by comparing the model predictions to experimental data on the static structure in section 3 and the dynamics of this system in section 4. The latter show a significant degree of cooperativity of the Li+ ions, which is characterized in section 5. Consequences for the existence of such a conduction mechanism in Li-ion conductor are drawn in section 6.

Vcharge =

∑ i,j>i

qiq j rij

(4)

In the case of silicate-based, phosphate-based, and sulfate-based materials, previous studies involving classical interaction potentials have either used formal charges for all the species49,50 or partial charges for the S and O atoms only.51,52 On the basis of the Bader atomic charges reported for LiMgSO4F by Ramzan et al.,47 which are close to the formal values, we have chosen the former approach. The dispersion potential includes two terms ⎡ C ij C ij ⎤ Vdispersion = − ∑ ⎢f 6ij (rij) 66 + f 8ij (rij) 88 ⎥ ⎢ rij ⎥⎦ rij i,j>i ⎣

(5)

where the short-range damping is described using the Tang− Toennies functions f ijn, which are of the form53 f nij

=1−e

−bDij rij

n



(bDij rij)k k!

k=0

(6)

The repulsion potential is modeled using a decaying exponential function Vrepulsion =

ij

∑ Aij e−B r

ij

(7)

i,j>i

A Gaussian term is also introduced to provide an extra steep repulsion at very short distances and therefore improve the stability of the potential.54−56 Finally, the induction term reads

2. METHODS 2.1. DFT Calculations. 2 × 2 × 2 supercells containing 16 LiMgSO4F units (128 atoms in total) were used for the DFT simulations, which were performed using the plane-wave DFT code CPMD.37 We employed the PBE exchange-correlation functional38 and a plane-wave cutoff of 1087.5 eV was found suitable to obtain sufficient convergence of the forces. We used pseudopotentials of the Troullier−Martins type for O, F, and S atoms39,40 and of the Goedecker−Teter−Hutter type for Li and Mg atoms.41,42 We have also determined the maximally localized Wannier functions43 for each configuration, which allowed us to calculate the induced dipole on each O and F atom.44 2.2. Interaction Potential Development. The accuracy of the results provided by molecular dynamics depends entirely on the reliability of the interaction potential employed. For some applications, it is possible to avoid this question by performing a DFT calculation at each time step (ab initio molecular dynamics45). Nevertheless, this approach is very expensive

Vinduction =

∑ [qiμαj g ij(rij) − q jμαi g ji(rij)]Tijα − ∑ μαi μβj Tijαβ i,j

+

i,j

∑ i

1 i2 |μ | 2α i

(8)

where α is the ion polarizability, and the Einstein summation convention is assumed. The set of induced dipoles {μ} are calculated at each time step by minimizing Vinduction. This differs from the approaches used in previous simulations of LiFePO4 and LiFeSO4, where the Drude oscillator method was used instead.49,52 A more noticeable difference is the inclusion of a short-range correction to the multipolar expansion,57 for which Tang−Toennies functions are also used: i

ij

4

g ij(rij) = 1 − c ij e−b rij ∑ k=0

18619

(bijrij)k k!

(9)

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C

Article

behavior for the Li+ ions over the time scale of our simulations; the material is probably not stable at such temperatures, but the use of classical MD (for which no chemical reaction can occur) renders possible such simulations. The time step for the integration of the equation of motions was set to 1.0 fs and the thermostat64 relaxation time to 20 ps. The total simulation time for each temperature was 8 ns. The electrostatics as well as dispersion interactions were calculated following the Ewald summation method.65−67 The calculations were performed on a 6 × 6 × 4 supercell containing 288 LiMgSO4F units. Bigger and smaller supercells were used to test the size dependence of the diffusion coefficients (see section 5). In the first series of simulations, we used the experimental lattice parameters of LiMgSO4F63 (cell 1 in the remaining text). It is worth noting that these parameters were obtained at room temperature. In order to analyze the effects of thermal expansion, which is likely to occur for the highest experimental temperature, we performed a second series of simulations with a larger cell. Here again, we have chosen to use the lattice parameters of LiFeSO4F since those correspond to a 4.5% larger volume than the ones of LiMgSO4F.6 The results obtained from the two sets of simulations (noted cell 1 and cell 2) will be discussed only when they differ in the remaining of the text.

These short-range effects for the dipoles have been shown to play an important role in ionic systems.48,58 In order to assign values to the various parameters introduced in the interaction potential, we use a generalized force-matching procedure59,60 where the reference data is obtained from firstprinciples DFT calculations on a series of ionic compounds (here, the fitting is performed on the forces applied on each atom and on the individual dipoles of O and F atoms). Here, the interaction potential cannot be considered as empirical since no experimental information is used in the parametrization procedure. This approach has also been shown to yield transferable potentials, which means that the parameters for a given pair interaction can be used in various environments.61 Some of the parameters were therefore kept equal to the value which had already been derived in other contexts (mainly for the Li−F and Mg−F pairs60). The others were determined by fitting the forces and dipoles obtained by DFT calculations on a series of typical LiMgSO4F configurations. As discussed by Wang and Van Voorhis, it is important to use structures originating from different sources for these reference data.62 Here, the DFT calculations were performed on two types of structures: on the one hand, we used configurations where the ions were positioned close to the sites determined in crystallographic experiments in a structural analogue,6 i.e., LiFeSO4F (the use of the experimental data63 of LiMgSO4F yields very small forces, which can hardly be fitted), and on the other hand, we used configurations that were generated from MD simulations (using parameters determined from a first fitting step). This yields a set of reference data spanning the important parts of the configuration space. An illustration of the quality of the fitted interaction potential is shown for the two typical structures in Figure 1. The agreement is

3. STRUCTURE OF LIMGSO4F The quality of the potential can now be assessed by comparing model predictions to experiments since no experimental data was used in the optimization of the model parameters. In this section, we will focus on static properties (lattice parameters, radial distribution functions, etc.). The building blocks of the LiMgSO4F structure are MgO4F2 octahedra and SO4 tetrahedra. Their arrangement consists of chains of corner-sharing octahedra, which are bridged by isolated tetrahedra,63 as shown in Figure 2. The lattice parameters of LiMgSO4 obtained from a geometry optimization are provided in Table 1, where they are also compared to experimental crystallographic data.63 An excellent agreement is found, with the largest differences of 0.022 Å for the cell parameter c and of 0.67° for the angle γ being well within 1.0% of their experimental counterpart. Such a level of agreement is usually obtained in DFT calculations, which illustrates the accuracy of the developed interaction potential. Next, in Figure 3, we report the radial distribution functions obtained from a short run at 300 K and with the experimental lattice parameters. The corresponding bond distances are reported in Table 2, together with the experimental data from Sebastian et al. The characteristic atomic distances also agree very well with experiments; in particular, we obtain the correct average S−O distance of 1.48 Å. The regions of maximum occupancy of Li+ ions, which are shown in blue in Figure 2, correspond well to the sites determined from X-ray diffraction data.63 In conclusion, we have validated the use of our interaction potential for studying the static properties of LiMgSO4. In the next section, we will focus on dynamic properties, such as Li-ion diffusion, which presents a greater challenge, as discussed in the introduction.

Figure 1. Illustration of the agreement between the components of the forces acting on the atoms of two LiMgSO4F configurations. Blue, obtained from the optimized interaction potential; red, obtained from DFT. Forces are given in atomic units. In configuration 1, the ions are placed on the crystallographic positions of a structual analogue (LiFeSO4F), while in configuration 2, their positions were generated from a short MD run.

overall excellent, which means that our potential should be of first-principles accuracy. A commented input file including all the fitted parameters is provided in the Supporting Information. 2.3. Molecular Dynamics Simulations. Molecular dynamics simulations were performed in the canonical (NVT) ensemble at nine temperatures ranging from 1100 to 1600 K. We investigated high temperatures in order to observe a diffusive

4. ELECTRICAL CONDUCTIVITY On the time scale of our simulations (8 ns), we observe a clearly diffusive behavior for the Li+ ions only for temperatures greater than 1100 K. For this reason, in the remainder of this article, we 18620

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C

Article

Figure 2. Snapshots of the LiMgSO4F structure along the x and y axes. SO4 tetrahedra and MgO4F2 octahedra are, respectively, shown in yellow and green. Regions of maximum occupancy of Li+ ions are shown in blue.

evolution of the positions of two Li+ ions at 1200 K is shown in Figure 4, where it appears clear that the diffusion occurs by

Table 1. Lattice Parameters with Experimental Values from Ref 63 a (Å) b (Å) c (Å) α (deg) β (deg) γ (deg)

simulation

experiment

5.1795 5.404 7.095 106.44 107.35 97.14

5.1623 5.388 7.073 106.68 107.40 97.81

Figure 4. Coordinates of two Li+ ions during a simulation performed at 1200 K .

successive hops. In agreement with previous static DFT17 and classical MD simulations34,52 performed on LiFeSO4F, we observe that a majority of these hops occur along the [111] channels of the tavorite structure; jumps from one channel to another occur less frequently. These hops appear to be strongly correlated (as shown by the concomitant change of the two particles in Figure 4), something that will be discussed later. A consequence of the jump mechanism is the existence of three régimes in the mean-squared displacements (MSD), calculated as discussed in the introduction and shown in Figure 5. For very short times, the Li+ ions follow a ballistic trajectory. Very rapidly, they bounce on the surrounding F and O atoms, and the MSD reaches a plateau value. The duration of this so-called caging régime, where the particles are temporarily trapped in cages formed by their neighbors, depends strongly on the temperature; typically, it lasts for around 100 ps at 1200 K (see Figure 5). For greater times, the ions start to hop on an adjacent site as depicted by the jumps in their coordinates in Figure 4. Consequently, the MSDs start to increase again, and their variation with time is linear. In their experimental study, Sebastian et al. reported measurements of the electrical conductivity in a wide range of temperatures.63 We remind that, in this material, there is no electronic contribution, so that the electrical conductivity equals the ionic conductivity. In molecular dynamics simulations, this quantity is normally defined as69

Figure 3. Partial radial distribution functions for selected cation−anion pairs.

Table 2. MD Bond Distances (Extracted from the Maximum of the First Peak of the Radial Distribution Functions in Figure 3) and Experimental Values Taken from Ref 63; Note That in the Case of the Li−O Distance, a Broad Range of Values Are Expected from Figure 3 ionic pair

MD

exptl

O−Mg F−Mg O−Li F−Li O−S

2.09 Å 1.95 Å 2.00 Å 1.82 Å 1.48 Å

2.05−2.12 Å 1.92−1.93 Å 2.00−2.19 Å 1.82−1.84 Å 1.46−1.48 Å

will focus on this high-temperature régime. It is worth noting that the internal energy varies linearly with temperature between 300 and 1200 K (see the Supporting Information), which corresponds to a constant heat capacity CV; therefore, there is no indication for a superionic transition68 in this system. The 18621

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C

Article

Figure 5. Mean-squared displacements (MSD) of the Li+ ions at 1200 K. The three components (X, Y, and Z) are also shown, and they are multiplied by a factor of 3 to facilitate the comparison.

σ=

βe 2 V

∫0

Figure 6. Individual and collective Li+ diffusion coefficients for the two different cell parameters.



J(t )dt

decreases with temperature; for example, in the first set of simulation conditions, it drops from 21.2 at 1100 K to 4.7 at 1600 K. Second, the cooperativy is more important in cell 1 compared to cell 2. Indeed, while in both cases, similar individual diffusion coefficients are obtained, there is a factor of approximately 2−3 between the two sets of data for the collective ones. The mechanism involved in this cooperative motion is further analyzed in the next section, which will allow us to interpret these results. As for the electrical conductivity, the behavior shown in Figure 6 indicates that the Nernst−Einstein approximation should not be used to determine it in LiMgSO4F. We have therefore calculated it using eq 11; the results for the two series of simulations are shown in Figure 7, where the experimental value

(10)

where β = 1/kBT, e is the electronic charge, and J is the charge current correlation function, but it is often more convenient to use the following equivalent expression:70,71 σ=

βe 2 1 lim ⟨| ∑ qiδ ri(t )|2 ⟩ V t →∞ 6t i

(11)

where δri(t) is the displacement of ion i, which carries a formal charge qi, at time t. The conductivity is also often approximated by the Nernst−Einstein expression σ NE = βe 2 ∑ ρi qiDi i

(12)

where the sum runs over all the mobile species of the system, with ρi being their number densities and Di the individual diffusion coefficient. In our case, only the Li+ ions diffuse, and this expression becomes

σ NE = βe 2ρLi DLi

(13)

By combining this expression with eq 2, we immediately see that σ and σNE are equivalent if ⟨δ ri(t ) ·δ rj(t )⟩ = 0, i ≠ j

(14)

i.e., if the correlations between the displacements of different ions are neglected. In our simulations, as observed in Figure 4, the hops of the Li+ ions within a [111] channel occur simultaneously, which means that this condition is not fulfilled. In order to check the validity of this expression, we have calculated the quantity D′Li =

1 1 lim ⟨| ∑ δ ri(t )|2 ⟩ N t →∞ 6t i

Figure 7. Comparison between the experimental (black triangles) and simulated (blue symbols) electrical conductivities for LiMgSO4F. The experimental points are extracted from ref 63.

from Sebastian et al.63 are also shown. Although the sets of data correspond to different temperature intervals, the extrapolation of the latter up to high-temperatures would give results in accord with our simulations.

(15)

This quantity can be seen as the collective diffusion coefficient, which is related to σ by the exact expression σ = βe 2ρLi D′Li

(16)

5. COOPERATIVE CONDUCTION MECHANISM Cooperative motion of atoms has already been observed in glassforming liquids for which a picture of dynamic facilitation has emerged.72 In such systems, it is common to detect the most mobile particles along the trajectory; in many examples, these appear to be spatially correlated.73,74 Here, following Vogel and Glotzer,75 we define a cluster of mobile Li+ ions as a group of the 10% most mobile particles (the ones with the largest mean-

The results obtained for the individual D and collective D′ diffusion coefficients of Li+ ions are summarized in Figure 6 for the two different cell geometries used in this study. The difference between them is spectacular since D′ is on average 1 order of magnitude larger than D. Such a result is the signature of a high degree of cooperativity in the Li+ ion diffusion. Two interesting features are also observed: first, the ratio D′/D 18622

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C

Article

squared displacements during a time interval Δt) that reside in the vicinity of each other. This means that a distance cutoff criterion has to be introduced for defining Li+ ions which are neighboring each other; on the basis of the typical lengths of hops as observed in Figure 4, we fixed this cutoff to 4.8 Å. The weightaveraged size of a cluster to which a mobile particle belongs is then given by S W (Δt ) =

∑n n2PS(n; Δt ) ∑n nPS(n; Δt )

(17)

where PS(n;Δt) is the probability distribution of finding a cluster of size n within a time interval Δt.75 Note that the use of a different cutoff value for the distance criterion leads to a change in the value of the maximum size of the cluster but not of the qualitative behavior. The results obtained for cell 1 are shown in Figure 8. Large cluster sizes are observed, which corroborates

Figure 9. Successive snapshots, taken along a [111] diffusion channel, highlighting the cooperative diffusion mechanism of Li+ ions at 1200 K. The SO42− tetrahedra are represented by the S−O bonds, while the Mg and F atoms are, respectively, shown as white and green spheres (note that the atoms are not shown using their usual van der Waals radii in order to enhance the Li+ ions positions). The Li+ ions appear as pink, red, blue, or purple spheres. The red Li+ ions jumps into the [111] channel, provoking a correlated displacement of all the blue ions, and finally the ejection of the purple one.

Figure 8. Weight-averaged cluster size, SW(Δt) for the mobile Li+ ions, at various temperatures (results obtained with the simulation cell 1).

with the high degree of cooperativity shown by the diffusion. As in glass-forming liquids, upon decreasing the temperature, the maximal cluster size increases, as well as the time interval for which it occurs. An example of collective diffusion behavior extracted from our simulation is shown in Figure 9. In this series of snapshots, the system is viewed along one of the [111] diffusion pathways. The SO42− tetrahedra are represented by the S−O bonds, while the Mg and F atoms are, respectively, shown as white and green spheres (note that the atoms are not shown using their usual van der Waals radii in order to enhance the Li+ ion positions). The Li+ ions appear as pink, red, blue, or purple spheres: the pink ones do not experience any jump during this portion of trajectory. The red one, which is in a different [111] channel at t = 240 ps, suddenly jumps, at t = 250 ps, in the channel where the blue Li+ atoms are. At that moment, the Li+ ions represented in blue start to propagate by hops from the left to the right; at each snapshot in the time interval between 250 and 280 ps, we observe a set of around 5 ions that are more closely packed, which shows the propagation of the excitation. The propagation stops at t = 290 ps, when the ion shown in purple is ejected in another [111] channel. In this example, if we consider the blue ions only, each of them has jumped by δrhop, while the overall charge was displaced by nine times δrhop, which highlights the difference between the individual and collective diffusion coefficients and the nonvalidity of the Nernst−Einstein relationship for systems with correlated motion.

We also note here that the cooperative character of the Li diffusion in this material has also important technical consequences. First, it appears clearly that the hypothesis that the Li+ ions are not interacting is not valid in this system, preventing the use of eq 1 (at least with the commonly used parameters). Second, special care must be taken when choosing the simulation box size. Indeed, the excitation shown in Figure 9 involves nine Li+ ions in the channel. This implies that simulation cells in which the channels contain at least twice this number of ions have to be used to avoid the spurious interaction of some Li+ ions with their periodic images. Indeed, when relatively small systems (96 LiMgSO4F units) were simulated, the calculated mean squared displacements were noisy, and the extraction of the diffusion coefficient was very difficult, even when trajectories of several tens of nanoseconds were used. We therefore performed calculations on system sizes ranging from 96 to 1600 LiMgSO4F units, corresponding to 768 to 12 000 atoms. We found that the collective diffusion coefficients converge as the system size is increased and are virtually system size independent when 6144 (768 LiMgSO4F units) atoms or more are simulated. Throughout this article, we have reported the results obtained from a supercell containing 288 LiMgSO4F units, corresponding to 2304 atoms (the 288 Li+ ions are shared between 12 channels of 24 ions each). The diffusion coefficients extracted from this 18623

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C



were very close to the converged values obtained from bigger systems, and this supercell allowed us to perform long simulation and thus to provide sufficient statistics at relatively low temperatures (1100−1300 K).

ASSOCIATED CONTENT

S Supporting Information *

Commented input file and plot of the internal energy vs temperature. This material is available free of charge via the Internet at http://pubs.acs.org.



REFERENCES

(1) Armand, M.; Tarascon, J.-M. Nature 2008, 451, 652−657. (2) Padhi, A. K.; Nanjundaswamy, K. S.; Goodenough, J. B. J. Electrochem. Soc. 1997, 144, 1188−1194. (3) Delmas, C.; Maccario, M.; Croguennec, L.; Le Cras, F.; Weill, F. Nat. Mater. 2008, 7, 665−671. (4) Sun, C.; Rajasekhara, S.; Goodenough, J. B.; Zhou, F. J. Am. Chem. Soc. 2011, 133, 2132−2135. (5) Nishimura, S.-I.; Nakamura, M.; Natsui, R.; Yamada, A. J. Am. Chem. Soc. 2010, 132, 13596−13597. (6) Recham, N.; Chotard, J.-N.; Dupont, L.; Delacourt, C.; Walker, W.; Armand, M.; Tarascon, J.-M. Nat. Mater. 2010, 9, 68−74. (7) Barpanda, P.; Chotard, J.-N.; Delacourt, C.; Reynaud, M.; Filinchuk, Y.; Armand, M.; Deschamps, M.; Tarascon, J.-M. Angew. Chem., Int. Ed. 2011, 50, 2526−2531. (8) Dambournet, D.; Chapman, K. W.; Chupas, P. J.; Gerald, R. E., II; Penin, N.; Labrugere, C.; Demourgues, A.; Tressaud, A.; Amine, K. J. Am. Chem. Soc. 2011, 133, 13240−13243. (9) Abouimrane, A.; Dambournet, D.; Chapman, K. W.; Chupas, P. J.; Weng, W.; Amine, K. J. Am. Chem. Soc. 2012, 134, 4505−4508. (10) Naille, S.; Jumas, J.-C.; Lippens, P.-E.; Olivier-Fourcade, J. J. Power Sources 2009, 189, 814−817. (11) Bhattacharyya, R.; Key, B.; Chen, H. L.; Best, A. S.; Hollenkamp, A. F.; Grey, C. P. Nat. Mater. 2010, 9, 504−510. (12) Key, B.; Morcrette, M.; Tarascon, J.-M.; Grey, C. P. J. Am. Chem. Soc. 2011, 133, 503−512. (13) Hung, I.; Pourpoint, F.; Grey, C. P.; Gan, Z. J. Am. Chem. Soc. 2012, 134, 1989−1901. (14) Ceder, G.; Chiang, Y. M.; Sadoway, D. R.; Aydinol, M. K.; Jang, Y. I.; Huang, B. Nature 1998, 392, 694−696. (15) Armstrong, A. R.; Lyness, C.; Panchmatia, P. M.; Islam, M. S.; Bruce, P. G. Nat. Mater. 2011, 10, 223−229. (16) Hautier, G.; Jain, A.; Ong, S. P.; Kang, B.; Moore, C.; Doe, R.; Ceder, G. Chem. Mater. 2011, 23, 3495−3508. (17) Mueller, T.; Hautier, G.; Jain, A.; Ceder, G. Chem. Mater. 2011, 23, 3854−3862. (18) Kim, J.; Middlemiss, D. S.; Chernova, N. A.; Zhu, B. Y. X.; Masquelier, C.; Grey, C. P. J. Am. Chem. Soc. 2010, 132, 16825−16840. (19) Castets, A.; Carlier, D.; Zhang, Y.; Boucher, F.; Marx, N.; Croguennec, L.; Menetrier, M. J. Phys. Chem. C 2011, 115, 16234− 16241. (20) Combelles, C.; Ben Yahia, M.; Pedesseau, L.; Doublet, M.-L. J. Phys. Chem. C 2010, 114, 9518−9527. (21) Barpanda, P.; Ati, M.; Melot, B. C.; Rousse, G.; Chotard, J.-N.; Doublet, M.-L.; Sougrati, M. T.; Corr, S. A.; Jumas, J.-C.; Tarascon, J.-M. Nat. Mater. 2011, 10, 772−779. (22) Van der Ven, A.; Ceder, G.; Asta, M.; Tepesch, P. D. Phys. Rev. B 2001, 64, 184307. (23) Morgan, D.; Van der Ven, A.; Ceder, G. Electrochem. Solid-State Lett. 2004, 7, A30−A32. (24) Amin, R.; Maier, J.; Balaya, P.; Chen, D. P.; Lin, C. T. Solid State Ionics 2008, 179, 1683−1687. (25) Maginn, E. J. Acc. Chem. Res. 2007, 40, 1200−1207. (26) Malik, R.; Burch, D.; Bazant, M.; Ceder, G. Nano Lett. 2010, 10, 4123−4127. (27) Marrocchelli, D.; Madden, P. A.; Norberg, S. T.; Hull, S. Chem. Mater. 2011, 23, 1365−1373. (28) Norberg, S. T.; Hull, S.; Ahmed, I.; Eriksson, S. G.; Marrocchelli, D.; Madden, P. A.; Li, P.; Irvine, J. T. S. Chem. Mater. 2011, 23, 1356− 1364. (29) Burbano, M.; Norberg, S. T.; Hull, S.; Eriksson, S. G.; Marrocchelli, D.; Madden, P. A.; Watson, G. W. Chem. Mater. 2012, 24, 222−229. (30) Chroneos, A.; Parfitt, D.; Kilner, J. A.; Grimes, R. W. J. Mater. Chem. 2010, 20, 266−270. (31) Kushima, A.; Parfitt, D.; Chroneos, A.; Yildiz, B.; Kilner, J. A.; Grimes, R. W. Phys. Chem. Chem. Phys. 2011, 13, 2242−2249. (32) Panchmatia, P. M.; Orera, A.; Rees, G. J.; Smith, M. E.; Hanna, J. V.; Slater, P. R.; Islam, M. S. Angew. Chem., Int. Ed. 2011, 50, 9328−9333.

6. CONCLUSIONS In this work, we studied the conduction mechanism of LiMgSO 4 F by means of extensive molecular dynamics simulations. First, we parametrized a dipole-polarizable interaction potential from first-principles DFT calculations, and then, we tested it by reproducing, to a high degree of accuracy, both structural and dynamical properties of LiMgSO4F. Our calculations were then used to characterize the conduction mechanism of this material. We showed that the diffusion of the Li+ ions occurs via correlated hops inside diffusion channels, which results in the much larger values calculated for their collective diffusion coefficient compared to the individual ones. These correlated motions are initiated by the jumps of ions from one channel to another. Although a three-dimensional diffusion is observed, the diffusion is therefore largely facilitated along the [111] channels. A consequence of the existence of such a cooperativity mechanism is that large simulation supercells are needed to correctly sample the whole diffusion process. Further work will focus on the application of the MD method to the study of the diffusion of Li+ ions in electronic conductors such as LiFePO4 and LiFeSO4F. This will require the introduction of additional features in the interaction potential. Indeed, the redox active species, Fe in those cases, should be allowed to change their redox state during the simulation. A promising way to introduce such effects in the framework of classical simulations is the approach where atoms from an electrode are held at a fixed potential. It has already been used to study the electrowetting of the metallic electrode surface by water76 or, more recently, to understand the charging mechanism of supercapacitors77 and could be extended to the study of Li-ion battery electrode materials.



Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (M.S.); D.Marrocchelli@ tcd.ie (D.M.). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS D.M. wishes to thank the Government of Ireland for an EMPOWER Postdoctoral Fellowship and the Trinity Centre for High Performance Computing for computing time. We wish to acknowledge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support (project tcche026b). G.W.W. would like to thank SFI for funding (grant numbers 08/RFP/MTR1044 and 09/RFP/MTR2274). M.S. thanks the theory group of the RS2E French network for fruitful discussions. 18624

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625

The Journal of Physical Chemistry C

Article

(33) Zhang, P.; Wu, Y.; Zhang, D.; Xu, Q.; Liu, J.; Ren, X.; Luo, Z.; Wang, M.; Hong, W. J. Phys. Chem. A 2008, 112, 5406−5410. (34) Adams, S.; Prasada Rao, R. Solid State Ionics 2011, 184, 57−61. (35) Yildirim, H.; Greeley, J. P.; Sankaranarayanan, S. K. R. S. Phys. Chem. Chem. Phys. 2012, 14, 4565−4576. (36) Lee, S.; Park, S. S. J. Phys. Chem. C 2012, 116, 6484−6489. (37) The CPMD Consortium. CPMD, version 3.13.2. http://www. cpmd.org. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (39) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 001993. (40) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 008861. (41) Goedecker, S.; Teter, M.; Hutter, J. Phys. Rev. B 1996, 54, 1703− 1710. (42) Krack, M. Theor. Chem. Acc. 2005, 114, 145. (43) Marzari, N.; Vanderbilt, D. Phys. Rev. B 1997, 56, 12847−12865. (44) Silvestrelli, P. L.; Parrinello, M. Phys. Rev. Lett. 1999, 82, 3308− 3311. (45) Vuilleumier, R. Lect. Notes Phys. 2006, 703, 223−285. (46) Yang, J.; Tse, J. S. J. Phys. Chem. A 2011, 115, 13045−13049. (47) Ramzan, M.; Lebègue, S.; Kang, T. W.; Ahuja, R. J. Phys. Chem. C 2011, 115, 2600−2603. (48) Madden, P. A.; Wilson, M. Chem. Soc. Rev. 1996, 25, 339−350. (49) Islam, M. S.; Driscoll, D. J.; Fisher, C. A. J.; Slater, P. R. Chem. Mater. 2005, 17, 5085−5092. (50) Armstrong, A. R.; Kuganathan, N.; Islam, M. S.; Bruce, P. G. J. Am. Chem. Soc. 2011, 133, 13031−13035. (51) Allan, N. L.; Rohl, A. L.; Gay, D. H.; Catlow, C. R. A.; Davey, R. J.; Mackrodt, W. C. Faraday Discuss. 1993, 95, 273−280. (52) Tripathi, R.; Gardiner, G. R.; Islam, M. S.; Nazar, L. F. Chem. Mater. 2011, 23, 2278−2284. (53) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726−3741. (54) Castiglione, M. J.; Wilson, M.; Madden, P. A. J. Phys.: Condens. Matter 1999, 11, 9009−9024. (55) Norberg, S. T.; Ahmed, I.; Hull, S.; Marrocchelli, D.; Madden, P. A. J. Phys.: Condens. Matter 2009, 21, 215401. (56) Burbano, M.; Marrocchelli, D.; Yildiz, B.; Tuller, H. L.; Norberg, S. T.; Hull, S.; Madden, P. A.; Watson, G. W. J. Phys.: Condens. Matter 2011, 23, 255402. (57) Stone, A. J. Theory of Intermolecular Forces; Oxford University Press: Oxford, U.K., 1996. (58) Wilson, M.; Madden, P. A. J. Phys.: Condens. Matter 1993, 5, 2687−2706. (59) Aguado, A.; Bernasconi, L.; Jahn, S.; Madden, P. A. Faraday Discuss. 2003, 124, 171−184. (60) Salanne, M.; Rotenberg, B.; Simon, C.; Jahn, S.; Vuilleumier, R.; Madden, P. A. Theor. Chem. Acc. 2012, 131, 1143. (61) Madden, P. A.; Heaton, R. J.; Aguado, A.; Jahn, S. J. Mol. Struct. 2006, 771, 9−18. (62) Wang, L.-P.; Van Voorhis, T. J. Chem. Phys. 2010, 133, 231101. (63) Sebastian, L.; Gopalakrishnan, J.; Piffard, Y. J. Mater. Chem. 2002, 12, 374−377. (64) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. J. Chem. Phys. 1992, 97, 2635−2643. (65) Karasawa, N.; Goddard, W., III. J. Phys. Chem. 1989, 93, 7320. (66) Aguado, A.; Madden, P. A. J. Chem. Phys. 2003, 119, 7471−7483. (67) Laino, T.; Hutter, J. J. Chem. Phys. 2008, 129, 074102. (68) Wilson, M.; Jahn, S.; Madden, P. A. J. Phys.: Condens. Matter 2004, 16, S2795−S2810. (69) Hansen, J.-P.; McDonald, I. Theory of Simple Liquids, 2nd ed.; Academic Press: New York, 1986. (70) Morgan, B.; Madden, P. A. J. Chem. Phys. 2004, 120, 1402−1413. (71) Salanne, M.; Simon, C.; Turq, P.; Madden, P. A. J. Phys. Chem. B 2007, 111, 4678−4684. (72) Keys, A. S.; Hedges, L. O.; Garrahan, J. P.; Glotzer, S. C.; Chandler, D. Phys. Rev. X 2011, 1, 021013. (73) Donati, C.; Douglas, J. F.; Kob, W.; Plimpton, S. J.; Poole, P. H.; Glotzer, S. C. Phys. Rev. Lett. 1998, 80, 2338−2341.

(74) Donati, C.; Glotzer, S. C.; Poole, P. H.; Kob, W.; Plimpton, S. J. Phys. Rev. E 1999, 60, 3107−3119. (75) Vogel, M.; Glotzer, S. Phys. Rev. Lett. 2004, 92, 255901. (76) Willard, A. P.; Reed, S. K.; Madden, P. A.; Chandler, D. Faraday Discuss. 2009, 141, 423−441. (77) Merlet, C.; Rotenberg, B.; Madden, P. A.; Taberna, P.-L.; Simon, P.; Gogotsi, Y.; Salanne, M. Nat. Mater. 2012, 11, 306−310.

18625

dx.doi.org/10.1021/jp304767d | J. Phys. Chem. C 2012, 116, 18618−18625