Effects of Li Confined Motion on NMR Quadrupolar Interactions: A Combined 7 Li NMR and DFT-MD Study of LiR2 (PO4 )3 (R=Ti or Zr) Phases.

Ab initio molecular dynamics (MD) simulations and NMR GIPAW (gauge including projector augmented wave) calculations have been used to analyze the coordination and mobility of Li ions in LiTi2 (PO4 )3 (rhombohedral), LiZr2 (PO4 )3 (triclinic), and LiZr2 (PO4 )3 (rhombohedral) phases. Significant discrepancies are observed between static calculations of 7 Li quadrupolar parameters and experimental values. The dynamical origin of this disagreement is demonstrated by incorporating in the calculations thermal vibrations and local motion of atoms with MD simulations. For LiTi2 (PO4 )3 , the quadrupolar constant associated with Li ions grows with temperature because the local symmetry of the system decreases, whereas for the Zr phases, the quadrupolar constant decreases because thermal vibrations reduce the anisotropy of the interaction. Finally, for both Zr phases, MD yields Li distributions that compare well with disorder reported from diffraction studies.

Introduction

LiR2(PO4)3 compounds with NASICON (Na SuperIonic CONductors) structure have been extensively studied as electrodes and electrolyte materials in Li batteries.1 Chemical stability and Li mobility are high in these compounds, making them good candidates for solid electrolytes in all‐solid‐state batteries (ASSB).2 Research into energy storage batteries has recently been deservedly recognized with the chemistry Nobel prize being awarded to A. Yoshino, M. S. Whittingham, and J. B. Goodenough.

The rhombohedral LiR2(PO4)3 compounds (R c space group) with NASICON structure are formed by R2(PO4)3 lanterns that share oxygen atoms with neighboring ones, building three‐dimensional networks.3 In such networks, M1 and M2 sites alternate along the conduction channels. In rhombohedral Ge or Ti phases, Li ions are sixfold coordinated at octahedral M1 sites4 (Figure 1 a). The incorporation of bigger R cations (Sn, Zr, or Hf) stabilizes crystal structures with triclinic C symmetry. In these phases, half of the cavities are occupied by two Li ions, whereas the other half of the cavities are empty. In triclinic structures, PO4 tilting favors the distorted fourfold coordination of Li ions at M12 sites near triangular bottlenecks that connect M1 and M2 cavities.5 In the triclinic structure of LiZr2(PO4)3 reported by Catti et al.,5c some disorder was detected, as the Li ions are located at two inequivalent sites with occupancies of 0.71 and 0.29 (Figure 1 c). Both sites have a symmetrical replica inside the same cavity.

Figure 1

(a) Ti and (b) Zr rhombohedral structures and (c) Zr triclinic structure. P orange, O red, Ti gray, Zr light blue, Li purple (white spheres partially colored in purple denote fractional occupancy).

Heating LiZr2(PO4)3 samples at 323 K causes a phase transition from triclinic to rhombohedral.5c, 6 In the rhombohedral NASICON structure, Li ions show disorder over six equivalent crystallographic sites with 1/6 occupancy (Figure 1 b). In this phase, Li ions are again located near triangular bottlenecks that connect M1 and M2 cavities.

NMR spectroscopy is often used to analyze cation sites. In the case of the 7Li magic angle spinning (MAS) NMR spectra, the quadrupolar constant (C Q) and the asymmetry parameter (η Q) allow the study of the distortion of the coordination polyhedron. The C Q experimental values obtained for different LiR2(PO4)3 phases and temperatures were previously reported7 (Figure 2). In LiZr2(PO4)3, C Q changes, as a consequence of the phase transition, from 186 kHz (triclinic structure) to 110 kHz (rhombohedral structure); whereas η Q changes from 0.3 to 0, respectively. In the Ti phase, an increase in C Q is observed when the temperature of the system is raised, whereas the asymmetry parameter remains equal to zero. The increment of C Q with temperature has been previously ascribed to the progressive occupation of M12 sites at the expenses of octahedral M1 sites.7

Figure 2

Variation of experimental quadrupolar constants of Ti and Zr phases with temperature.

The chemical shift anisotropy (CSA) and the electric field gradient (EFG) are physical magnitudes that can be represented by second rank tensors. Such tensors can be visualized as ellipsoids whose axes are defined by the principal values of the associated 3×3 matrixes. The quadrupolar constant and the asymmetry parameter are related to the EFG tensor through the following expressions [Eq. (1)]:

where V, V, and V are the principal values of the EFG tensor at the nucleus being studied and Q is its quadrupole moment. Although calculated C Q values can be positive or negative, experiments only provide absolute values, therefore, comparisons in this paper will be done through absolute values.

Discrepancies between calculated values and experimental measurements of chemical shift anisotropies, dipolar interactions, and quadrupolar interactions have often been assigned to vibrations of atoms induced by temperature. Dračínský and Hodgkinson8 reported the effects of dynamics on the calculated CSA and EFG tensors for a set of amino acids and nucleic acid bases. Changes in the orientation of these tensors owing to atom vibrations result in more spherical average tensors, which in turn give smaller chemical shift anisotropy and quadrupolar constant values than those derived from static model tensors. Fluctuations of less than 10° on the principal axes caused a decrease of 2–5 % in the absolute values of the quadrupolar coupling constants.8 Other sources of disagreement can be poor crystallographic models and/or disorder in synthesized samples.

Conduction pathways and long‐range Li/Na mobility in NASICON compounds have been studied by using the nudged elastic band method.9 Ab initio10 and classical11 molecular dynamics (MD) methods also have been applied in a few cases to study Li diffusion in NASICON electrolytes with different compositions to the ones considered in this work.

In the present work, DFT calculations have been performed to obtain optimized structures for Ti and Zr phases. As the calculated quadrupolar 7Li MAS NMR parameters for these static structures differ considerably from the experimental values, we used MD to incorporate dynamical effects owing to Li mobility into the calculations, which restores the agreement between theory and experiment.

Results and Discussion

Owing to the narrow 7Li chemical shift range and the similarity between structures of the studied compounds, local motions are not expected to cause significant variations in the isotropic chemical shift values. However, the quadrupolar interaction is very sensitive to the distribution of charges around the targeted nucleus and it is expected to be affected by dynamical effects. Consequently, we will focus this work on quadrupolar interaction. Static 7Li spectra and fittings of LiTi2(PO4)3 and LiZr2(PO4)3 phases at high and low temperatures are shown in Figure 3.

Figure 3

Static 7Li spectra (blue) and simulations (red) of LiTi2(PO4)3 phase at 300 and 450 K and LT and HT LiZr2(PO4)3 phases.

The accuracy achieved in the simulations discards the possibility of important disorder or secondary phases in the samples. Experimental values for the quadrupolar constants and the asymmetry parameters are summarized in Table 1.

Table 1

Experimental and theoretical 7Li quadrupolar (C Q and η Q) constants for the different phases of this study.

	Phase/model	7Li |C Q| [kHz]	7Li η Q
Experimental	Ti (300 K)	42	0
	Ti (450 K)	61	0
Calculated	optTiAatiq	21	0
	optTi303	19	0
	optTi448	21	0
Experimental	Zr (300 K)	186	0.3
Calculated	optZrCatti99	249	0.34
Experimental	Zr (425 K)	110	0
Calculated	optZrCatti00	154–221	>0.85
	optZrM1	27	0

Ground state structures

Our starting points are the optimized geometries for the Ti phase, the low‐temperature (LT) Zr phase, and the high‐temperature (HT) Zr phase.

Ti phase

The crystallographic structure of the rhombohedral LiTi2(PO4)3 phase reported by Aatiq et al.12 (which showed atomic forces up to 0.61 eV Å−1) was geometry optimized, that is, positions of the atoms were iteratively improved until the maximum force on any of them was lower than an imposed threshold of 0.01 eV Å−1. The unit cell parameters were fixed to experimental values (structure optTiAatiq). NMR GIPAW (gauge including projector augmented wave) calculations were performed for this structure. Quadrupolar parameters obtained for 7Li are C Q=21 kHz and η Q=0 (Table 1). The calculated asymmetry parameter agrees with the experimental value. However, the calculated quadrupolar constant is smaller than the experimental value. Extrapolating experimental data to T=0 K improves the agreement, which already indicates the importance of taking into account the temperature (Figure 2). To study the dependence of quadrupolar parameters with the size of the unit cell, the Ti phase has been optimized with parameters obtained in our lab (for details, see the Supporting Information, Section S2) at 303 and 448 K (structures optTi303 and optTi448, respectively). The maximum atomic forces of these structures prior to the optimization procedure were similar to the previous one. In the structure optTi303, where the unit cell is slightly smaller than that reported by Aatiq, we have obtained a smaller quadrupolar constant than for the model optTiAatiq (19 kHz vs. 21 kHz). Hence, further contractions could explain the experimental value (17 kHz) at 135 K (smallest value in Figure 2). The structure optTi448 has almost the same unit cell parameters as the one reported by Aatiq and we have obtained similar quadrupolar parameters. An optimization of the structure proposed by Aatiq, allowing unit‐cell parameters to vary, undergoes an expansion of 3.9 % in cell volume. The calculated quadrupolar constant of this model is 35.5 kHz. The preferential elongation of the cavity in the z direction is probably the responsible for this increase in the C Q value. Although a unit cell expansion induced by temperature can be proposed as a cause for higher quadrupolar constants, it cannot be considered as the only reason for the increment observed in experimental values. Inclusion of thermal effects will be further studied below.

The distance between the Li ions and their six nearest O atoms is 2.10 Å for the R c model displaying the smallest cell parameters (optTi303) and 2.13 Å for the structures optTiAatiq and optTi448. These values are close to the sum of the ionic radii of Li and O atoms (2.16 Å), which justifies the sixfold coordination for lithium ions in this phase.

LT Zr phase

In the structural model reported by Catti et al.5c for LiZr2(PO4)3 at room temperature, Li ions are located at two sites with occupancies of 0.71 and 0.29. To carry out the geometry optimization of this structure, only the Li site with the highest occupancy was considered to host the Li ion. In the resulting optimized ground state structure (optZrCatti99), Li ions are shifted 0.085 Å from the reported site.5c The maximum atomic forces diminish from 5.3 eV Å−1 in the model without optimization to <0.01 eV Å−1 after optimization. Li distances to the four nearest O atoms were in the range 1.98–2.27 Å, suggesting a planar fourfold coordination. The best fit plane of the four coordinated O atoms has a maximum deviation of 0.15 Å and the distance from the Li ion to this plane is 0.25 Å. 7Li quadrupolar parameters obtained for the optimized LT Zr phase (optZrCatti99) are summarized in Table 1. The calculated C Q value is −249 kHz, its absolute value is 34 % larger than the experimental one (186 kHz), whereas the calculated asymmetry parameter agrees with the experimental one. The second largest axis of the EFG ellipsoid points to the nearest O atom, which is located at less than 2 Å from the Li atom, and its largest axis points to a point between the other Li atom of the cavity and the further Zr atom.

HT Zr phase

In the LiZr2(PO4)3 rhombohedral phase reported by Catti and Stramare,13 the Li−O distances would be 2.50 Å if the Li ions were located at M1 sites (Figure 4 b). In contrast, if Li ions move towards any of the six triangular windows related by −3 axes, their coordination number will lower and three Li−O distances will become shorter. In the reported structure,13 Li ions are disordered over six equivalent positions near the triangular windows that connect M1 and M2 cavities (Figure 1 b). The Li−O distances to the three nearest O atoms of the cavity are 1.92, 2.43, and 2.45 Å. An O atom outside the cavity is located at 2.29 Å, suggesting a four‐fold coordination. To perform the structural optimization, the symmetry of the system was reduced from R c to P1 and one out of six Li sites in each cavity was chosen to host the Li ion (Figure 4 a). As our model has two cavities and, depending on the choice of Li site to be occupied, six inequivalent models were built. These models have maximum atom forces in the range 1.68–2.79 eV Å−1. After the geometry optimization, the maximum atom forces are <0.01 eV Å−1 and the final energy of the considered models (ensemble optZrCatti00) remains in a narrow range of 0.03 eV f.u.−1. Small differences detected in the energy calculations indicate that, although the position of Li ions can be affected by Li ions in the neighbor cavities, all proposed Li configurations are accessible at room temperature. The optimized structures locate all Li ions at 2.0–2.1 Å from the three O atoms forming the nearest window of the M1 cavity.

Figure 4

HT Zr models: (a) Li atom near the centroid of a triangular window and (b) Li atom at M1 site. In both cases, the calculated EFG tensors are rendered in blue. The scale factor for the EFG tensor in structure (b) is twice the scale factor in structure (a).

7Li quadrupolar constants obtained for the optimized HT Zr phase (optZrCatti00) range between 154 and 221 kHz with asymmetry parameters higher than 0.85 (Figure 4 a). Comparing these results with the experimental values (C Q=110 kHz and η Q=0), there are significant discrepancies for both C Q and η Q. Therefore, an alternative model with Li ions located at the center of M1 cavities was considered (Figure 4 b). In this case, the geometry optimization yields a structure (optZrM1) with an energy only 0.06 eV f.u.−1 higher than the most stable configuration, although it displays long Li−O distances (2.39 Å). The 7Li quadrupolar constant obtained for this model is 27 kHz. This value confirms the high symmetry of M1 sites, but still significantly differs from the experimental value. Calculated EFG ellipsoids corresponding to Li ions located at the triangular windows that surround M1 and to Li ions located at M1 sites are depicted in Figure 4 a and b, respectively. For the structure with the Li ion located near a triangular window (Figure 4 a) the EFG ellipsoid with disc‐like shape is almost coplanar with the plane defined by the Li atom and the two nearest Zr atoms and one of its two largest principal axis almost points to the nearest Zr atom. For the structure with the Li ion located at the M1 position, the largest axis of the ellipsoid is aligned with both Zr atoms. The Zr atoms seem to be crucial for the shape and orientation of the 7Li EFG ellipsoids.

As Li ions mostly display local mobility inside the M1 cavity, we can create new five symmetry equivalent structures for optZrCatti00. Averaging the EFG tensors for the Li atoms in these structures, we can estimate the effect on the quadrupolar interaction of Li ions jumping from one window to another. Following such a procedure, the value obtained for the quadrupolar constant is 171.5 kHz and for the asymmetry parameter 0. Therefore, we conclude that the calculated value for the quadrupolar constant is too large to be successfully compared with the experimental one. The proposed averaging is likely wrong because it implies that all the conformations are equally probable and thermal vibrations of atoms are ignored.

Molecular dynamics

To check the effect of thermal motions on the 7Li NMR parameters, ab initio molecular dynamics simulations were performed (see movies in the Supporting Information). For all the studied phases, the P−O bonds display lengths with a Gaussian type distribution centered at 1.54–1.55 Å and with a full width at half height of 0.08–0.09 Å. Zr−O and Ti−O bonds display length distributions centered at 2.06 and 1.92 Å, respectively, with full widths at half height of 0.16–0.18 Å. However, Li−O bonds display a more complex behavior, which is explained in detail in the next paragraphs for each phase.

Two trajectories of 5 ps were calculated at 300 K and two others at 450 K, taking as starting models the optimized structures optTi303 and optTi448. Movies showing one trajectory at each temperature can be found in the Supporting Information. Lithium ions have a larger mobility and explore larger volumes than other atoms. The shape of such volume is an oblate ellipsoid, with its minor axis oriented along the −3 axis direction of the NASICON structure. Figure 5 a and b show the superposition of 250 frames taken every 20 fs from trajectories calculated at 300 K and at 450 K. The disorder observed at 450 K is clearly higher than at 300 K.

Figure 5

Superposition of frames (one cavity) taken from the trajectories every 20×10−15 s (Li: purple spheres, O: red, Ti: gray, and Zr: light blue) for the Ti phase at 300 K and 450 K (a and b) and for the HT Zr phase at 425 K (c). Occupancy density of Li ions at a certain distance from: (d) the centroid of the six O atoms delimiting the cavity and (e) those six O atoms.

Figure 5 d shows the occupancy density of Li ions at increasing distances from the centroid of the six O atoms forming the M1 cavity. The distribution of distances is wider for 450 K than for 300 K, and its maximum moves further away from the center of the cavity with temperature (see the Supporting Information, Section S2 for further discussion). At 300 K, Figure 5 d shows that the probability of finding Li ions at a maximum distance of 0.43 Å is 90 %, whereas at 450 K only 70 % probability is within that range. The maximum of occupancy density also displaces: at 300 K it is 0.25 Å, whereas at 450 K it is 0.33 Å. At 450 K, the separation between the most distant Li positions is 1.9 Å. Figure 5 e shows that distances between the Li ions and the six nearest O atoms are in the range between 1.75 and 2.90 Å for 300 K and in the range between 1.70 and 3.50 Å for 450 K (the most frequent distances being in the range 2.1–2.3 Å in both cases).

Single snapshots of any of the studied trajectories usually show a maximum atomic deviation of 0.5 Å from ideal R c symmetry positions. This deviation is mainly related to Li ions. The maximum atomic deviation of the framework in each frame of the trajectory is usually larger than 0.2 Å at 300 K. Calculating the average positions of atoms for different simulation times we find that, for 1 ps, the maximum deviation of the structure stays always below 0.1 Å and if we do not consider Li ions, it stays below 0.03 Å. The average positions for a 5 ps simulation decrease the maximum deviation to 0.02 and 0.01 Å considering the whole structure or the framework. We can even postulate smaller deviations for longer simulation times. As deviations from symmetry can be embodied in thermal factors, this result is consistent with the rhombohedral model obtained by diffraction studies. Deviations are bigger for the 450 K trajectory, than for the 300 K trajectory. In the high‐temperature simulation, we find maximum deviations for the whole structure as large as 0.13 Å and 0.05 Å for 1 ps and 5 ps averaging times, respectively.

Two trajectories of 5 ps were calculated at 300 K, taking optZrCatti99 as the starting model. A movie showing one of these trajectories can be found in the Supporting Information. Details of the Li coordination and occupancy density of Li ions versus Li–O distance in this phase are given in the Supporting Information, Section S3.

Figure 6 shows the superposition of 250 frames taken every 20 fs from a 300 K trajectory of the LT Zr phase. The clouds of the Li ions are elongated and located near the centroids of opposite triangular windows that connect M1 and M2 cavities. These positions are near the crystallographic site displaying the highest occupancy. Most of the Li ions can be found within a 0.7 Å sphere centered in that site. The cloud of Li ions is aligned with the direction connecting the two disordered positions found by diffraction techniques, although the dispersion of Li positions predicted by MD does not reach the second site reported in the diffraction model. The most distant positions occupied by the same Li ion during this simulation are 1.6 Å away from each other. The average positions for a 5 ps simulation show a maximum deviation of 0.02 Å from ideal positions displaying the C symmetry found by diffraction methods.

Figure 6

Superposition of 250 frames (every 20 fs) of one LT Zr phase MD trajectory at 300 K (Zr: blue, O: red, and Li: purple). The two positions detected by powder neutron diffraction are represented by white spheres at the centers of two gray spheres of 0.7 Å diameter.

Three trajectories of 10 ps were calculated at 425 K. In two MD simulations, the starting model was the one with the lowest energy of the optZrCatti00 ensemble of structures. In the third calculation, the starting model was the one with the second lowest energy. A movie showing one of these trajectories can be found in the Supporting Information. No net difference was found between the three trajectories. All of them show the system transforming between the proposed conformations.

Figures 5 c and 7 show the superposition of 500 frames from one of the 10 ps trajectories. Li ions are represented as purple spheres spreading from the center of the M1 cavity and beyond the windows limiting the cavity. Although Li ions concentrate in six symmetry equivalent directions, the cloud does not display R c symmetry. By comparing Figure 5 b and c, we can see that Li ions in the HT Zr phase explore a wider volume than in the Ti phase at comparable temperatures. Figure 7 shows the six O atoms (red) that limit the M1 cavity and the two Zr atoms (light blue) attached to them at their mean positions in a 10 ps trajectory. The small green spheres (numbered from 1 to 6) are the geometric centroids of each triangular window. Measuring the distance between the Li ion and these six centroids, the residence time of the mobile ion near each window can be estimated.

Figure 7

Detail of one M1 cavity. Superposition of 500 frames (every 20 fs) of one HT Zr phase MD simulation at 425 K. Zr: blue, O: red, and Li: purple (Zr and O are located at mean positions of the trajectory). Small green spheres numbered from 1 to 6 are located at the centroids of the triangular windows that connect M1 and M2 cavities.

Figure 8 represents the distance from one Li ion to the nearest window centroid at each moment in one of the trajectories. After the equilibration time, the Li ion is located near window 6. After 372 fs, it moves closer to window 5 for a short time and afterwards, it jumps to window 4. After several further jumps, this trajectory finishes at t=10 ps with the Li ion near window 3.

Figure 8

Li hopping between equivalent triangular windows of the M1 cavity is represented by vertical dashed lines. Li vibrations at each window are also depicted.

Analyzing the evolution of all Li ions involved in the three trajectories, we find that the mean residence time at one window is 190 fs (53 jumps, on average, occur in 10 ps). Most of the jumps take place between neighboring windows, but if the Li ion approaches the centroid of the cavity, the jump can occur between more distant windows. For all the trajectories studied, the six windows around M1 sites were visited by Li ions, but the sum of residence times in each of them is not homogeneous, it differs by as much as a factor two in the simulated intervals. For the trajectory shown in Figure 8, the Li ion spends as much as 268 fs in window 2 and only 107 fs in window 6. We expect that for longer simulation times, these differences will be attenuated, as the probability of finding the Li ion near any of them is equivalent. Owing to the computational effort needed to perform this type of ab initio simulations, we have not calculated longer trajectories, which would decrease the dispersion of calculated values, but should not significantly affect the averages themselves, as convergence of average C Q values is already obtained for 10 ps trajectories (see convergence study in Section S4 in the Supporting Information). The occupancy density of Li ions at a certain distance of the six nearest O atoms is represented in Figure 5 e. The sharp maximum for the three nearest ones is at 2.05 Å, whereas for the rest of the O atoms the graph shows a wider distribution of distances centered at 3.1 Å. The most distant positions occupied by the same Li ion during this simulation are 4.4 Å away from each other. Figure 5 d shows that Li ions in the Ti phase remain close to the centroid of the M1 cavity, whereas for the HT Zr phase the probability of finding the Li ions at this centroid is scarce. Its wide maximum is centered 1 Å away from this site. The best way to describe the cloud of Li ions with discrete Li sites would be the use of six symmetry‐equivalent positions near the windows connecting M1 and M2 cavities as is done in the reported diffraction model.

We have investigated atom deviations from nominal R c positions in one of the trajectories of the HT Zr phase. Structures generated by averaging atom positions after 1 ps and 10 ps simulation times give deviations from 0.5 to 0.33 Å considering Li ions and deviations from 0.18 to 0.04 Å if Li ions are not included. These values are clearly higher than the ones found for the Ti phase at comparable simulation times. In the Ti phase, averaging over 1 ps provides deviations from R c symmetry comparable to the ones found for the HT Zr phase for 10 ps. The total simulated span calculated in this study is not enough to obtain an average structure that complies well with the symmetry obtained by diffraction in this phase, but it describes quite well the Li confined motion in the M1 cavity of the NASICON structure.

Summarizing the results obtained by MD simulations, we observe a tendency for Li to coordinate with O atoms with distances shorter than 2.3 Å. For small M1 cavities (Ti phase) this means positions near the centroid of the cavity, whereas for bigger cavities (HT Zr phase) this means positions near the centroids of triangular windows. No Li jumps between cavities were observed at the temperature and simulation times used in this work. However, simulations obtained are useful to evaluate the effect of dynamics on the quadrupolar interaction, as we prove in the next section.

Averaging 7Li quadrupolar interactions

If the residence time of a given nucleus at different structural sites is much shorter than the time required for NMR detection (1/ω 0), the observed quadrupolar interaction will correspond to the averaged electrical field gradient felt by that nucleus during the acquisition time (ms range). As motions observed in previous simulations are fast with respect to the NMR timescale, the prediction of 7Li quadrupolar C Q and η Q constants requires the averaging of the electric field gradient felt by the Li ions through the calculated trajectories. As a word of caution, we notice that our dynamics simulations only cover a few ps, whereas experimental NMR averaging involves longer times. Therefore, we assume a fair sampling is obtained in the ps range, and the system will merely be revisiting similar configurations at longer times.

GIPAW calculations were performed on structures taken from the four simulated trajectories: two at 300 K and two at 450 K. These structures were extracted every 20 fs, resulting in 250 structures for each 5 ps trajectory. Following this procedure, 250 EFG tensors for four independent Li atoms (A and B from trajectory 1 and C and D from trajectory 2) were determined at each temperature. C Q calculated values go from 12 to 224 kHz at 300 K and from 14 to 342 kHz at 450 K along the trajectories. To consider the effect of atom vibrations on quadrupolar parameters, the nine components of the EFG tensors obtained for each inequivalent Li ion were averaged. The four resulting tensors were diagonalized and quadrupolar parameters calculated. Values obtained for each temperature are summarized in Table 2. These ellipsoids have their major axes oriented along the −3 axis direction of the NASICON structure, in a similar way to that shown in Figure 4 b.

Table 2

C Q and η Q values obtained after averaging 250 EFG tensors calculated for snapshots of trajectories for the Ti phase at 300 and 450 K.

	300 K			Symm. av.			450 K			Symm. av.
	|C Q|	η Q		|C Q|	η Q		|C Q|	η Q		|C Q|	η Q
Li(A)	47	0.20		47	0.00		71	0.15		71	0.00
Li(B)	48	0.42		46	0.00		62	0.11		61	0.00
Li(C)	46	0.10		45	0.00		67	0.33		66	0.00
Li(D)	47	0.13		46	0.00		65	0.23		64	0.00
Mean value	47	0.21		46	0.00		66	0.19		66	0.00
Experimental	42	0.0					61	0.0

Li(A) and Li(B): first trajectory; Li(C) and Li (D): second trajectory. Columns labeled symm. av. show results obtained after symmetry averaging. Numbers in bold compare experimental and calculated average values.

Calculated quadrupolar constants considering MD are higher than those obtained from static ground state structures (47 vs. 19 and 66 vs. 21 for 300 and 450 K, respectively). It is shown that vibrations of the framework and motion of Li ions in M1 cavities have the global effect of lowering the local symmetry around Li ions. The tendency of higher C Q at higher temperatures is confirmed by previously reported experiments on this phase.14 This increment is related to the expansion of the cell, which promotes Li occupancy of less symmetrical sites as the volume of the cavity becomes larger and to a larger thermal motion of lithium ions. In Table S5 in the Supporting Information, the effects of cell expansion and thermal motion at 300 and 450 K are considered separately. Calculated quadrupolar constants considering MD are only 8–12 % higher than the experimental ones and confirm the dependence of C Q with temperature. However, the short averaging times used (5 ps) do not perfectly reproduce the axial symmetry observed in rhombohedral phases (η Q=0). To take into account dynamical effects over longer times, the random movement of Li ions within the cavity and the ternary axis of the structure, five equivalent ellipsoids were generated considering a −3 axis and, a final averaging of the six ellipsoids was performed. These calculations further reduce C Q values by 0–2 kHz and the asymmetry parameters become equal to 0 (Table 2, symm. av.).

EFG tensors were similarly averaged for the two trajectories calculated for the Zr phase at 300 K; values obtained are summarized in Table 3. In this phase, framework vibrations and the motion of Li ions changes the value of the quadrupolar constant from −249 to −204 kHz. The average absolute value is 10 % higher than the experimental one (186 kHz). The effect of the temperature in this phase is a reduction in the anisotropy of the interaction, which lowers the value of the quadrupolar constant. This result is opposite to the previous case, where Li ions are located at high symmetry sites and thermal motion lowers the symmetry increasing the C Q.

Table 3

C Q and η Q values obtained by averaging 250 EFG tensors calculated for snapshots of trajectories for LT Zr phase at 300 K.

		|C Q|	η Q
	Li(A)	200	0.37
	Li(B)	213	0.32
	Li(C)	204	0.36
	Li(D)	197	0.34
	Mean value	204	0.35
	Experimental	186	0.3

Li(A) and Li(B): first trajectory; Li(C) and Li (D): second trajectory. Numbers in bold compare experimental and calculated average values.

In this case, the motion of Li ions consists of vibrations around the centroids of triangular windows that connect M1 and M2 cavities and hopping between the six equivalent windows that surround M1 sites. To perform the average, 500 frames were extracted from each trajectory (one every 20 fs). Values calculated for C Q on these snapshots cover a wide range: going from a couple of tens to 400 kHz. The averaged quadrupolar parameters obtained for each one of the three 10 ps trajectories are given in Table 4.

Table 4

C Q and η Q values obtained by averaging 500 EFG tensors calculated for snapshots of trajectories for the HT Zr phase at 425 K.

	425 K		Symm. av.
	|C Q|	η Q	|C Q|	η Q
Li(A)	136	0.14	136	0.00
Li(B)	143	0.18	141	0.00
Li(C)	130	0.17	130	0.00
Li(D)	141	0.09	140	0.00
Li(E)	132	0.17	132	0.00
Li(F)	135	0.20	135	0.00
Mean value	136	0.16	136	0.00
Experimental	110	0.0

Li(A) and Li(B): first trajectory; Li(C) and Li(D): second trajectory; Li(E) and Li(F): third trajectory. Columns labeled symm. av. show results obtained after symmetry averaging. Numbers in bold compare experimental and calculated average values.

The C Q value is lowered from values in the range 168–221 kHz, obtained for the static configurations of the system, to 136 kHz. This value is 24 % higher than the experimental value. If we visualize the orientation of the averaged ellipsoids (see figure in the table of contents for one Li atom of this phase), we observe that their longer axes are almost aligned with the ternary axis of the R c NASICON structure. Averaging the EFG tensors lowers the asymmetry parameter from values higher than 0.85 to values lower than 0.2. This result suggests that the simulation time should be further increased to well reproduce quadrupolar interactions with axial symmetry (η Q=0). As commented above, the cloud of Li ions does not display R c symmetry. To override the lack of R c symmetry over short simulations times, we introduce five symmetry‐equivalent trajectories to those calculated. Thus, five symmetric EFG tensors equivalent to those previously calculated were generated and a final averaging of each group of six tensors was performed (column labeled symm. av. in Table 4). Such a procedure gives C Q values similar to those previously calculated, suggesting that the simulation describes well the effect of motions that lower the anisotropy of the quadrupolar interaction.

Finally, we carried out MD simulations for the rhombohedral phases at the same temperatures considered above, but with structural models obtained after geometry optimization of reported structures allowing unit cell parameters to vary. The optimized structures have a tendency to expand the volume of the cell to minimize stresses related to the details of the simulation (e.g., the choice for the exchange‐correlation functional). For the Ti phase, MD simulations show Li ions exploring larger regions of space, partly because M1 cavities are larger. Therefore, the quadrupolar constants obtained by averaging the calculated parameters in snapshots of these trajectories are larger than those reported in previous paragraphs. In the case of the HT Zr phase, the time spent in each window increases and the number of jumps between windows decreases. The region near the M1 site has a lower probability of being visited by the Li ions. The quadrupolar constants obtained by averaging the results calculated on snapshots along the new trajectories are lower than those previously reported for the HT Zr phase. The volume of the cavity seems to be crucial for defining the volume explored by Li ions and thus the estimation of quadrupolar interaction. Arbi et al.15 reported an increment of the quadrupolar constant as the unit cell augments along the LiTi2−Zr(PO4)3 series. This experimental observation is another manifestation of the relation between M1 cavity volume and quadrupolar interaction in NASICON compounds.

The high crystallinity of the studied phases facilitates setting up theoretical models. In addition, the Li atoms have labile coordination and they are located within a cavity, which they are unlikely to leave. These two facts make them ideal systems for analyzing the confined motion of Li ions and their effects on quadrupolar coupling. Studies of this type could be carried out on other well characterized materials in which there is confined movement (in the ps range) of an atom. If the system had a moderate degree of positional/occupational disorder, it could be modeled with supercells or a set of simple models.16 However, the computing effort would be higher than for ordered models. Ab initio molecular dynamics, which has served in this study to analyze the confined mobility of lithium ions, could be used to study the conduction paths of ions in materials of interest in lithium or post‐lithium batteries.

Conclusions

DFT‐NMR calculations on optimized structures significantly disagree with experimental values for 7Li NMR quadrupolar parameters. Ab initio molecular dynamics simulations of LiTi2(PO4)3 and LiZr2(PO4)3 show short‐range Li mobility on the ps scale. For the Ti phase, Li ions move near M1 sites, whereas for the LT Zr phase they move near M12 sites and for the HT Zr phase they exchange between M12 sites within the same cavity. Simulations tell us that Li ions tend to be coordinated to O atoms at distances shorter than 2.3 Å. Trajectories calculated for HT and LT Zr phases justify the Li disorder found in neutron diffraction studies: the clouds of Li ions can be represented by a small number of partially occupied Li sites with large thermal factors.

Our results bring agreement between theory and experiment and provide a physical explanation for the apparent disagreement with static calculations: averaging dynamical effects are crucial to get correct values for C Q and η Q. Atom vibrations and Li mobility cause a decrease in the 7Li C Q in phases where Li ions are located at low symmetry sites (HT and LT Zr phases). On the contrary, MD simulations predict an increment of 7Li C Q with temperature in the Ti phase, where Li ions occupy a high symmetry site. Comparing static and dynamical calculations, we find a decrease in C Q of 18 % (LT Zr phase). Such a difference is larger than any reported in previous works analyzing dynamical effects on quadrupolar parameters.8, 17 Moreover, in the systems previously studied no increment in the C Q was found, whereas we reproduce for the Ti phase the experimentally observed increment in C Q of more than 100 % at 300 K and more than 200 % at 450 K. In the case of the HT Zr phase, MD simulations not only predict a considerable decrease in C Q (from 19 % to 38 % depending on the starting model used in the MD simulation), but simultaneously explain a drastic change in the η Q value: from about 1 for static structures to about 0 owing to the dynamics of the system. Although long‐range Li motion above 330 K in the Ti phase has been reported,14 our simulations show that at the considered temperatures the main dynamical effects on the 7Li quadrupolar interaction are caused by the confined thermal motion of Li ions, which is more ubiquitous and likely than jumps between cavities.

Experimental Section

Synthesis and characterization

LiTi2(PO4)3 and LiZr2(PO4)3 compounds were prepared by following the reported method.5c, 14 X‐ray diffraction (XRD) was used to assess purity of the prepared compositions. XRD patterns were recorded with the CuKR radiation (λ=1.5405981 Å) in the temperature range 303–698 K, by using a high‐temperature AP HTK10 camera adapted to a PW 1050/25 Phillips diffractometer. Powder XRD patterns were collected in the 10–70° 2θ range with a step size of 0.02° and a counting time of 0.1 s/step. The Fullprof program18 was used to determine the unit cell parameters of prepared samples. 7Li MAS NMR spectra were recorded with an AVANCE 400 Bruker spectrometer (9.4 T) at 155.50 MHz. Spectra were obtained after π/2 pulse irradiation (3.5 μs) with a recycling time of 15 s. The number of scans was 8–400. The fitting of the NMR spectra was carried out with the Dmfit software.19 The fitting of spectra was performed by trial and error procedures. Mean errors on C Q and η parameters were below 3 kHz and 0.1, respectively.

Theoretical formalism and simulations

Calculations were carried out by using Density Functional Theory (DFT)20 as implemented in CASTEP21 (Materials Studio Suite 2017R2). The formalism employs the gauge including projector augmented wave (GIPAW)22 algorithm, enabling the reconstruction of the all‐electron wave function in the presence of a magnetic field. Core–valence interactions were described by using ultrasoft pseudopotentials23 and the generalized gradient approximation (GGA) PBE24 functional was used. Geometry optimizations were performed with a kinetic energy cutoff of 1200 eV and a Monkhorst–Pack25 k‐point grid spacing of 0.04 Å−1. Optimization of atomic coordinates were pursued until the energy difference, maximum atomic force, and maximum atomic displacement fell below 5×10−6 eV/atom, 1×10−2 eV Å−1, and 5×10−4 Å tolerances, respectively. For GIPAW quadrupolar parameter calculations, integrals over the Brillouin zone were performed by using a k‐point grid spacing of 0.07 Å−1. Wavefunctions were expanded in planewaves with kinetic energy up to a cutoff energy of 1400 eV, which was found to be enough to get converged C Q values. Li mobility was investigated with Born–Oppenheimer molecular dynamics simulations (BOMD)26 in the NVE ensemble.27 The microcanonical is the appropriate ensemble for a problem where the number of particles and the volume is fixed and at the same time it provides a simple way to keep track of the effect of rounding errors while integrating the equations of motion by monitoring the conservation of the energy. CASTEP21 has been used to perform these simulations, with the following parameters: an integration time step of 1×10−15 s, “on the fly” pseudopotentials, an energy cutoff for plane waves of 600 eV, and a Monkhorst–Pack grid with a k‐point spacing of 0.07 Å−1. The extended Lagrangian formulation was used.28 Two trajectories of 5×10−12 s were calculated for the Ti phase at 300 and 450 K, two more trajectories of 5×10−12 s were calculated for the LT Zr phase and three trajectories of 10×10−12 s were calculated for the HT Zr phase. The initial 200 fs, involved in the equilibration of the system, were disregarded. To perform the averaging of quadrupolar interactions, snapshots were extracted every 20×10−15 s. GIPAW NMR calculations were performed on every selected snapshot and each component of the EFG tensor was averaged to obtain the average tensor. From this average tensor, the quadrupolar constant and the asymmetry parameter were obtained by matrix diagonalization, according to the expressions [Eq. (1)].

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.

Supplementary

Click here for additional data file.

Supplementary

Click here for additional data file.

Supplementary

Click here for additional data file.

Supplementary

Click here for additional data file.