INVESTIGATION INTO THE DYNAMIC AND STRUCTURAL PROPERTIES OF THE LITHIUM GARNET SERIES ( Li 7 - x La 3 Zr 2 - x Ta x O 12 , x = 0 - 2): A COMBINED MOLECULAR DYNAMICS AND QUASI ELASTIC NEUTRON SCATTERING STUDY By Matthew Klenk A DISSERTATION Submitted to Michigan State University i n partial fulfillment of the requirements f or the degree of Materials Science and Engineering Doctor of Philosophy 2019 AB STRAC T INVESTIGATION INTO THE DYNAMIC AND STRUCTURAL PROPERTIES OF THE LITHIUM GARNET SERIES ( Li 7 - x La 3 Zr 2 - x Ta x O 12 , x = 0 - 2): A COMBINED MOLECULAR DYNAMICS AND QUASI ELASTIC NEUTRON SCATTERING STUDY By Matthew Klenk The lithium garnet series Li 7 - x La 3 Zr 2 - x Ta x O 12 (x = 0 2) has shown great promise as a solid electrolyte material, however the room temperature conductivity is currently too low to find wide commercial success. In order to better understand the mechanisms of ionic diffusion within the crystal, a combined molec ular dynamics and quasi - elastic neutron scattering (QENS) study was investigated. Using molecular dynamics simulations, we are able to easily probe atomic scale events that are usually difficult to examine experimentally. Like the local arrangement of lith ium on its sublattice, or how the trajectory of lithium ions is affected by the nearest neighbor capable of capturing both the residence time and mean jump distance experimentally allowing us to directly compare experimental and simulated intermediate scattering functions I,(Q,t). Overall, we saw good agreement between the two techniques, both predicting a jump - diffusion model in the form described by Singwi and Sjöla nder. Three different simulation models were employed in this study, two using classical molecular dynamics (MD), while a third using density functional theory (DFT) based calculations. All three model types are used to first better understand the phase tr ansformation behavior for the end member composition Li 7 La 3 Zr 2 O 12 , which undergoes a characteristic phase transformation from an ordered tetragonal to disordered cubic phase at 900 K. First DFT methods are used to better understand what role the selection of an electron exchange - correlation functional plays on the accuracy of lattice parameter and phase transformation behavior. In total 14 different functional forms are investigated. Similarly, two different classical MD models, one - el, where each atomic nucleus is connected by a spring potential to an - model that treats each atom as a point charge, which can be used for larger and faster s imulations. The dynamics of the two end member compositions Li 5 La 3 Ta 2 O 12 (L5LT) and Li 7 La 3 Zr 2 O 12 (L7LZ), were looked at using the core - shell model with respect to properties like the conductivity, self - diffusivity, thermodynamic correction factor, and ent ropy of configuration. While the core - only model is used to investigate the finite - size effects of atomic simulation, by changing the number of particles within the simulation by using four different crystal sizes for L7LZ. Simulation cells consisting of 1 ×1×1 (192 atoms), 2×2×2 (1536 atoms), 3×3×3 (5184 atoms), and 4×4×4 unit cells (12288 atoms) were generated in order to find convergent behavior to the properties highlighted above. Having determined a 3 x 3 x 3 simulation provides adequate accuracy, a ver ity of garnet compositions corresponding to [Li] = 5, 5.5, 6, 6.25, 6.5, 6.75, & 7 were generated to determine the optimal composition for use as an electrolyte material. Our simulations predict that the best performing room temperature composition corresp onds to when [Li] = 6.5 corresponding to the maximum lithium concentration that results in a disordered cubic phase at room temperature. Lastly, we look at the role lithium disorder plays in the phase transformation behavior of L7LZ, and the use of excess entropy calculations as a means of determining the performance of an electrolyte material. iv Dedicated to Friends and family v ACKNOWLEDGEMENTS I would first like to thank my advisor Dr. Lai for the wonderful opportunity granted to me to peruse this work at Michigan State University . You have been a wonderful mentor to me, showing me the hard work and dedication required to be a scienti st in today research environment. I have been in constant awe of the way he carries himself professionally and personally , always providing a platform for discussion, and a willingness to be open to new ideas of exploration. If I am so lucky to follow a similar path , I could not ask for a better person to provide a model for my behavior . Your kindness , understanding, and willingness to accept the difficulties I endured during this process has forever changed me and I could not ask for a better Ph D advisor. Thank you from the bottom of my heart, I cou ld not be where I am today without your help. I would like to thank my colleagues , former and current, for the all the times we have shared together over the years. Without the guidance and work of Yuxing Wang, I be gin to know where to start in the work presented in this thesis . You are instrumental in my success and I will always be grateful for our discussions that lead to many of the ideas explored here . Rengaranjan Shanmugam, you ar e one of the kindest people I have been so lucky to meet. I f I ever felt drained or ran out of ideas, talking with you always seemed to bring a fresh perspective and reinvigorated me to dive back in and solve any issue troubling me. Junchao Li, Jin Dai , an d Q ian Chen , I wish you all the success in the world with your coming endeavors. Jin yo ur fresh perspective on our work continually opened me up to new possibilities and areas of study , and your assistance in performing tasks with me on the garnet series has undoubting led to some of the results presented here. Chen and Junchao , your work on other materials systems has vi kept me captivated and motived to continue working in materials science . I would like the thank the undergraduate res earch assistants Sydney Boeberitz, Dong Feng, Christina Lopez, and Yalun Cai . Each one of you contributed directly to the work presented in this thesis and without your assistance, it would not be possible. Lastly, I would like to thank all my friends and family for all of your support throughout the years. Mom and Dad, y ou never doubted me for second , and gave me the confidence to follow my dreams. Without your love and support, I would not be the person I am today . My sisters Mary and Becky, thanks for putting up with my ramblings on the state of my research . Having you there for me whenever I needed it cleared my head and always reminded me of unbelievable girlfriend Kasey. Y ou always see the potential in me and give me the structure and motivation to be my best self . The neutron scattering experiments used resources at the Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory. I wish to acknowledge the Michigan State University High Performance Computing Center and the Institute for Cyber - Enabled Research for access to their computing resources. The work is financially supported by the Ceramics Program of National Science Foundation (DMR - 1206356) and partially by the James Dyson Foundation F e llowship . vii TABLE OF CONTENTS L IST OF TABLES ................................ ................................ ................................ .......................... x LIST OF FIGURES ................................ ................................ ................................ ....................... xi CHAPTER 1: Introduction ................................ ................................ ................................ ............. 1 1.1 Motivation for study ................................ ................................ ................................ .............. 1 1.2 The lithium garnet crystal ................................ ................................ ................................ ..... 8 1.3 Tetragonal to cubic transformation in Li 7 La 3 Zr 2 O 12 ................................ ........................... 12 1.3.1 Experimental observations ................................ ................................ ........................... 12 1.3.2 Simulations of garnet in literature ................................ ................................ ................ 13 1.4 Experimental measurement of conductivity and diffusivity ................................ ............... 16 1.4.1 Conductivity ................................ ................................ ................................ ................. 16 1.4.2 Diffusivity ................................ ................................ ................................ ..................... 17 CHAPTER 2: Materials and m ethods ................................ ................................ ........................... 19 2.1 Sample preparation ................................ ................................ ................................ .............. 19 2.2 Quasi - elastic neutron scattering ................................ ................................ .......................... 19 2.3 Molecular dynamic simulations ................................ ................................ .......................... 20 2.3.1 Core - shell model ................................ ................................ ................................ .......... 21 2.3.2 Core - only model ................................ ................................ ................................ ........... 23 2.4 Finite size effects in L7LZ ................................ ................................ ................................ .. 25 2.5 Effect of Ta/Zr substitution ................................ ................................ ................................ . 26 2.6 Density functional theory: XC functional analysis ................................ ............................. 26 2.6.1 Background ................................ ................................ ................................ ................... 26 2.6.2 Simulation method ................................ ................................ ................................ ........ 30 CHAPTER 3: Molecular dynamics analysis methods ................................ ................................ .. 32 3.1 Diffusivity and conductivity calculation methods ................................ .............................. 32 3.2 Calculation of collective diffusivity ................................ ................................ .................... 34 3.3 Practical considerations for fast - ion conductors ................................ ................................ . 37 3.4 Dynam ................................ ....................... 39 3.5 Thermodynamic factor ................................ ................................ ................................ ........ 47 3.6 Excess entropy ................................ ................................ ................................ ..................... 52 3.7 Statement on structure of this thesis ................................ ................................ .................... 55 CHAPTER 4: Effects of electron exchange - correlation functionals on the density functional theory simulation of phase transformation of fa st - ion conductors: A case study in the Li garnet oxide Li 7 La 3 Zr 2 O 12 ................................ ................................ ................................ ........................ 58 4.1 Abstract ................................ ................................ ................................ ............................... 58 4.2. Results and discussion ................................ ................................ ................................ ........ 59 4.2.1. Phase behavior at 1000 K: all 14 functionals ................................ .............................. 59 4.2.2. Lattice parameters at all temperatures: 8 functionals ................................ .................. 61 viii 4 .2.3 Correlation between XC functionals and phase behavior/volume ............................... 66 4.3 Conclusion ................................ ................................ ................................ ........................... 67 CHAPTER 5: Core - shell Modeling of L5LT and L7LZ ................................ .............................. 68 5.1 Abstract ................................ ................................ ................................ ............................... 68 5.2 Results and discussion ................................ ................................ ................................ ......... 68 5.2.1 Validation of simulation process ................................ ................................ .................. 68 5.2.2 L7LZ van Hove correlation functions (VHCF) ................................ ............................ 74 5.2.3 Lithium distribution and its relation to phase transition ................................ ............... 78 5.2.4 Lithium clusters ................................ ................................ ................................ ............ 82 5.3 Conclusions ................................ ................................ ................................ ......................... 91 CHAPTER 6: Lithium self - diffusion in a model lithium garnet oxide Li 5 La 3 Ta 2 O 12 : a combined quasi - elastic neutron scattering and molecular dynamics study ................................ ................... 92 6.1 Abstract ................................ ................................ ................................ ............................... 92 6.2 Results and discussion ................................ ................................ ................................ ......... 93 6.2.1 Static structure factor ................................ ................................ ................................ .... 93 6.2.2 Overall intermediate scat tering function and dynamic structure factor ....................... 93 6.2.3 Species - resolved dynamics: La, Ta, and O ................................ ................................ .. 97 6.2.4 Li dynamics ................................ ................................ ................................ .................. 97 6.3 Conclusions ................................ ................................ ................................ ....................... 102 CHAPTER 7: Finite - size effects on the molecul ar dynamics Simulation of fast - ion conductors: A case study of lithium garnet oxide Li 7 La 3 Zr 2 O 12 ................................ ................................ ........ 10 3 7.1 Abstract ................................ ................................ ................................ ............................. 103 7.1.1 Background ................................ ................................ ................................ ................. 103 7.2 Results and discussion ................................ ................................ ................................ ....... 104 7.2.1 Phase characterization ................................ ................................ ................................ 104 7.2.2 Self - diffusivity, ionic conductivity, and Haven ratio ................................ ................. 107 7.2.3 Thermodynamic factor an d Fickian diffusivity ................................ .......................... 110 7.3 Conclusions ................................ ................................ ................................ ....................... 113 CHAPTER 8: Modeling Li 7 - xLa 3 Zr 2 - x Ta x O 12 using core only model ................................ ........ 114 8.1 Abstract ................................ ................................ ................................ ............................. 114 8.2 Results and discussion ................................ ................................ ................................ ....... 115 8.2.1 Phase characterization of 3x3x3 LLZT cells ................................ .............................. 115 8.2.2 Ionic conductivity, ................................ ................................ ................................ ...... 117 8.2.3 Self - diffusivity ................................ ................................ ................................ ............ 119 8.2.4 Haven ratio, thermodynamic factor ................................ ................................ ............ 120 8.2.5 Li site occupancy and regular solution model ................................ ............................ 124 8.2.6 Li sub - lattice and configurational energy ................................ ................................ ... 127 8.3 Conclusions ................................ ................................ ................................ ....................... 132 CHAPTER 9: Future work and conclusions ................................ ................................ ............... 134 9.1 Future work ................................ ................................ ................................ ....................... 134 9.1.1 Exploration in to multivalent doping schemes ................................ ........................... 134 ix 9.1.2 Electrolyte - electrode interface modeling. ................................ ................................ .. 135 9.1.3 Evaluation of density correlation functions ................................ ................................ 135 9.2 Conclusions ................................ ................................ ................................ ....................... 136 BIBLIOGRAPHY ................................ ................................ ................................ ....................... 140 x L IST OF TABLES Table 1: Core - shell Buckingham and spring pair interaction parameters for LLZT .................... 22 Table 2: Core - only force field parameters ................................ ................................ .................... 25 Table 3: Exchange and correlation enhancement factors ................................ ............................. 29 Table 4: Properties derived b y changing the weighting factor in Equation 16 ............................. 36 Table 5: Coherent and incoherent scattering powers of elements involved in the presen t study. 46 Table 6: Summary of simulation results using different XC functionals ................................ ..... 60 Table 7: Parameters for the regular solution model ................................ ................................ .... 125 Table 8 : Refined regular solution parameters ................................ ................................ ............. 127 xi LIST OF FIGURES Figure 1: Diagram of a simplified electrochemical cell in the c harged state. The ion (green) is transported through the electrolyte, the electron (yellow) travels through an external circuit, where it recombines with the ion in the cathode. ................................ ................................ ............ 2 Figure 2: Conductivity of various electrolyte materials compared to that of LLZT 2 . The blue box designates the room temperature conductivity for the various materials ................................ ....... 6 Figure 3 : Space filling model of L5LT (left), La (green) and Ta (brown) make a framework structure for lithium (blue) to diffuse. (right) Lithium tetrahedral (pink) and octahedral (blue) sites in the lithium garnet (Wang 35 ). ................................ ................................ ............................... 8 Figure 4 : Diagrams of lithium occupancy for three different compositions of LLZT at 300 K. (a) Li 3 La 3 Te 2 O 12 , (b) Li 5 La 3 Ta 2 O 12 , and (c) Li 7 La 3 Zr 2 O 12 . ................................ ............................... 10 Figure 5: Schematic showing the relationship between possible lithium Wyckoff positions in the cubic (left) and tetragonal (right) crystals. Filled circles are occupied lithium sites and hollow circles are unoccupied sites. ................................ ................................ ................................ .......... 11 Figure 6: Mean squared displacement of lithium ions in L7LZ at various temperatures during a 1000 ps simulation. ................................ ................................ ................................ ....................... 33 Figure 7: The Li Li velocity autocorrelation i n the longitudinal (Black) and transverse (red) scattering directions. ................................ ................................ ................................ ..................... 37 Figure 8 : QENS HWHM fitting of L5LT at 600 K for possible diffusion models. ..................... 43 Figure 9: Finite size e ffects of - 1 for L7LZ at 1400 K for 1 x 1 x 1 to 4 x 4 x 4 simulation cells. (a) Absolute magnitude of - 1 as a function of the inverse length of the sub - volume, solid lines denote guides to the eye. (b) - 1 as a function of the volume ratio of the sub - volu me and the simulation cell with a third order polynomial fit. ................................ ................................ ......... 48 Figure 10: Simulated Li Li P DF for L7LZ at 1400 K for various simulation cell sized. ........... 51 Figure 11: Time dependent lattice parameter for the 14 tested functionals at 1000 K. ................ 59 Figure 12: Time dependent lattice param eters for 8 XC functionals at 900 K ............................. 62 Figure 13: Time dependent lattice parameter for 8 XC functionals at 300 K ............................... 63 Figure 14: Time dependent lattice parameters for 8 XC functionals at 700 K ............................. 63 Figure 15: Average lattice parameters for 8 XC functionals at 300 K, 700 K, 900K, and 1000 K (red). Experimental lattice parameter by X - ray diffraction by Matsuda. ................................ ..... 65 xii Figure 16: Plot of exchange enhancements factors of GGA functionals employed in the study. own in the dotted line. ........................... 66 Figure 17 : (a) Lattice parameter of L5LT compared to neutron and x - ray diffraction studies. (b) Atomic pair distribution function from total scattering experiments compared to MD simulations ................................ ................................ ................................ ................................ ....................... 70 Figure 18: (a) conductivity for L5LT and (b) L7LZ as a function of temperature. Blue lines represent the conversion of MSD to D*, and then applying the Nernst - Einste in (N - E) equation (Equation 8). Red lines are the thermodynamically corrected ( - 1*N - E) Nernst - Einstein equation. (c) Comparison between the inverse thermodynamic factors ( - 1) for L5LT and L7LZ using the fluctuation method (Equation 38). ................................ ................................ ................. 72 Figure 19: (a) Self - part of the van hove correlation function and (b) its spatial and temporal projection at 300 K. (c) Distinct - part of the van Hove correlation function and (d) its spatial and temporal projection at 300 K. ................................ ................................ ................................ ....... 74 Figure 20: (a) Mean squared displacement (MSD) for LLZ at 300 K and 1100 K. (b) Partial pair distribution function of Li Li pairs at different temperatures compared to that of Adams et al. 75 Figure 21: (a) Self - part of the van Hove correlation function and (b) its spatial and temporal projection at 1100 K. (c) Self - par t of the van Hove correlation function and (d) its temporal projection of a continuous random - walk model at 1100 K. ................................ .......................... 76 Figure 22: Lithium density maps on (111) for LLZ. (a) 300 K by Revitfield refinement, (b) 300 K, (c) 500 K, (d) 800 K , and (e) 1100 K by MD simulation. Isosurface level is 0.3 Å - 3 ............. 78 Figure 23: (a) Tetrahedral and octahedral cage occupancy for L7LZ as a function of temperature. Error bars represent the degree of dynamic fluctuation in the occupancy observed during MD simulation. Values from Bernstein et al are shown for comparison. (b) MSD for selected atoms occ upying three different cages type at 300 and 500 K. ................................ ............................... 79 Figure 24: Lithium density maps for L5L T derived from MD simulation at (a) 700 K, (b) 475 K, and (c) 300 K along the [100] direction. Isosurface level of 0.1 Å - 3.Two - dimensional Li density maps for the (111) plane, with a distance of 23 Å to the origin, cutting through Td - Oh - Td cages at (d) 700 K, (e) 475 K, and (f) 300 K with isosurface levels from 0 to 1 Å - 3. ........................... 81 Figure 25: Possible l ocal arrangements of Li with respect to their nearest neighbors ................ 82 Figure 26 : Average number of each type of possible lithium nearest neighbor arrangement as a function of temperature for L5LT centered upon Td (a) and Oh (b) sites, LLZ for Td (c) and Oh (d). ` ................................ ................................ ................................ ................................ ............... 83 Figure 27: Observed dynamical events at 500 K for L5LT. Yellow circles represent Li ions. Pink squares are tetrahedral sites and blue rectangles are octahedra lc ages. In each example, the local environment, Li trajectories, and geometries of bottlenecks as Li goes through the faces are illustrated. ................................ ................................ ................................ ................................ ...... 85 xiii Figure 28: Examples of lithium dynamics at 300 K for (a) T14 cluster, (b) a T04 cluster, and at 1100 K for a (c) T14 cluster, and (d) a T04 cluster. Yellow dots represent Li atoms. Pink squares and blue rect angles schematically represent Td and Oh cages ................................ ..................... 87 Figure 29: Histograms of lithium position o n the face of the bottle neck with respect to distance from the nearest oxygen or shared edge of oxygen polyhedral for LLT and LLZ. (a) diagram of bottle neck with solid lines representing distance to nearest oxygen, dashed lines representing distance to n earest face. Shown are events for LLT at (b) 550 K, (d) 900 K, and (f) 1100 K. Shown for LLZ at (c) 600 K, (e) 900 K, and (g) 1100 K. ................................ ............................. 90 Figure 30: (a) The static structure factor obtained from the QENS experiment and MD simulations at 700 K. The results are vertically offset for clarity and lines serve as a guide to the eye. (b) Neutron - weighte d (including both coherent and incoherent intermediate scattering functions calculated from MD at 700 K for several Q ................................ ................................ .. 93 Figure 31: (a) Dynamic structure factor S(Q, E) from QENS experiments at 700 K for different Q. (b) Dynamic structure factor S(Q,E) from QENS at Q = 0.45 Å - 1, for different temperatures. (c) An example of the Wens fit with background, delta, and Lorentzian functions convoluted with the resolution function for 700 K and Q = 0.35 Å - 1, along with the residuals of the fit. ( d) HWHM ( ) of the Lorentzian as a function of Q2 at different temperatures. Solid lines are the fits to the Singwi - Sj öl ander model. ................................ ................................ .............................. 95 Figure 32: (a) Intermediate scattering functions of species La, Ta, and O at 700 K from Md for - Waller formula. (c) Comparison of the mean square displa cement Uiso from MD with that obtained from neutron powder diffraction. 3 ................................ ................................ ................................ ................................ ... 96 Figure 33 : (a) intermediate scattering functions of Li at 700 K from MD for selected Q. The solid lines are the fit of the KWW model to data after the first 10 ps. (b) Q - dependence of the KWW stretching parameter at different temperatures. (c) Schematic of the near est - neighbor tetrahedral (pink) and octahedral (blue) cages in lithium garnet oxides. Filled yellow circles are the ................................ ................................ 98 Figure 34: (a) HWHM ( ) experimentally measured using QENS and equivalent of Li diffusion derived from MD simulations in L5LT at 700 K. Solid lines are the fit of the SS model . (b) Self - diffusivity of Li in lithium garnet oxides of different compositions at different temperatures using various measurement probes. The activation energies for both MD and experimental sources are also provided. ................................ ................................ ....................... 99 Figure 35: (a) Mean residence time ( ) of Li diffusion in Li5La3Ta2O12 at sites obtained from MD simulation and QENS experiments. Solid lines are an Arrhenius fit. The residence time obtained from SLR - NMR experiments are also shown for Al - doped Li7La3Zr2O12 (cubic, C - L7LZ). (b) Mean jump leng th of Li diffusing in Li5La3Ta2O12 obtained from MD simulation and QENS experiments. Solid lines are guides to the eye. (c) Schematic of the relation between crystallographic sites for Li within tetrahedral (Td) and octahedral (Oh) cages in Li5La3Ta2O12 sh owing key distances from previous results. ................................ ................................ ............. 100 xiv Figure 36: Finite - size effects on the phase tra nsition. (a) Average lattice parameters as a function of temperature for 1×1×1 4×4×4 simulations. Solid lines are guide to the eye. Experimental data from Matsuda et al 98 and Larraz et al 43 are shown for comparison. (b) Fluctuation of lattice paramet ers as a function of number of atoms for selected temperatures. Solid lines are the linear fit (slopes close to - 0.5). (c) Time evolution of lattice parameters for different cell sizes at 900 K. ................................ ................................ ................................ ................................ ................. 105 Figure 37: Change in lattice parameter for 3x3x3 and 4x4x4 under heating and cooling cycles 106 Figure 38: Self - diffusivity as a function of temperature for four different simulation sizes. Literature data from the FPMD simulation (1×1×1) by Jalem et al are shown for comparison. The inset shows details at high temperatures. Solid lines are guide to the eye. ........................ 107 Figure 39: Fi nite - size effects on the ionic conductivity for 3×3×3 and 4×4×4 cells. Values from experimental measurements by Wang et al 1 and Matsui et al4, along with FPMD s ................................ ...... 109 Figure 40: Finite - size effects on the Haven ratio for four simulation cell sizes. Solid lines are guide to the eye. ................................ ................................ ................................ .......................... 110 Figure 41: Finite - - 1 at different temperatures. Solid lines are guide to the eye. ................................ ................................ ................................ .. 111 Figure 42: Finite - size effects on the Fickian diffusivity for four simulation sizes. .................... 112 Figure 43: Anisotropic NPT simulations on 3x3x3 simulation cells. (A) Unit cell lattice parameters for LLZT series. Solid l ines are guides to the eye. (B) Unit cell volume at 300 K, 600 K, 900 K, and 1200 K for LLZT series. Experimental data from Logeat et al, Mukhopadhyay et al, and Wang et al. are shown for comparison. (C) Time dependent lattice parameter for L6.75LZT at 300 K, 400 K, 500 K, 600 K. ................................ ................................ ................ 115 Figure 44: (a) Arrhenius plot of conductivity for LLZT series. Squares are calculated conductivity values; lines are best fit from 1000 K 500 K then extrapolated down to 300 K. Experimental total conductivty for LLT shown by Thangadurai. Inset shows high temperatur e conductivity with lines as guides. (b) Iso - temperature plot of conductivity as a function of lithium concentration. 1200 K, 900 K, 600 K are direct measurements, 300 K is extrapolated from best fit line in (a). (c) Activation energy as a function of lith ium concentration for LLZT series Compared to reported activation energies for various garnet compositions with standard deviation bars. ................................ ................................ ................................ ............................. 117 Figure 45: (a) Arrhenius plot of self - diffusivity for LLZT series with extrapolated best fit lines to 300 K. Inset shows high temperature 1400 K to 1000K with lines acting as guides to the eye. (b) Isothermal plot of self - diffusivity as a function of lithium concentration. 300 K data taken from extrapolated lines of (a). ................................ ................................ ................................ ............. 119 Figure 46: (a) Haven ratio for Li diffusion in the LLZT series at 300 K, 600 K, 900 K, and 1200 K and (b) for Li Li pairs at 300 K, 600 K, 900 K, and 1400 K . ............................ 120 xv Figure 47: for (a) Li - La pairs, (b) Li - Ta pairs, (c) Li - O pairs, and (d) Li - Zr pairs. ................................ ................................ ................................ ................................ ............ 123 Figure 48: (a) tetrahedral and (b) octahedral occupancy as a function of Li concentration and temperature. ................................ ................................ ................................ ................................ 125 Figure 49: Normalized excess entropy for the LLZT series as a function of composition and tempe rature. ................................ ................................ ................................ ................................ 129 Figure 50: Lithium nuclear density maps for LLZT series along the [100] direction a t 300 K, 600 K, 900 K, and 1200 K. Isosurface level 0.3 Å - 3 . ................................ ................................ ......... 131 1 C HAPTER 1 : Introduction 1.1 Motivation for s tudy O ne of the more transformative technologies of the past 30 years has been the lithium - ion battery 6 - 10 . When first introduced to power wha t we would now consider rudimentary phones, it would have been inconceivable to predict the wide range of applications where we can now find the technology. The untethering of electrical devices from outlets has become a defining characteristic of modern s ociety , the ability to use computers, smartphones, or tablets anywhere has led to fundamental shift s in the ways we work and play . The ubiquity of mobile computing and communication has even made possible to connect to every other human on the planet by as 11 The technology that underpins and powers the se innovation s is the lithium - ion battery, and the development of new battery technology will help usher in a new age of environmentally friendlier and safer person al devices . The lithium - ion battery revaluation goes beyond just communication devices or computers. We can now not only take electrical power with us anywhere, the push for the development of affordable electric vehicles will propel us to that anywhere. 12 - 13 This major shift in fundamental transportation technology may help pave a path to a greener future by reducing the use of fossil fuels and help in curbing global warming. Major developments in battery capacity and safety will be needed to achieve such lofty goals but these are challenges being heavily researched today. Because lithium - ion batteries are so widely used it is paramount that the technol ogy provides a safe interface with the user to ensure best experience possible. Newspaper and television articles have been filled with stories like the exploding smartphone or that of the 2 es work without fault. In the cases above, a flammable liquid organic mixture is used transport lithium - ions between two electrodes. W hen this liquid is overcharged, heated, or decomposed through aging , catastrophic consequences can occur in the form of th ermal runaway. 14 - 15 A s eries of irreversible exothermic reactions that can cause a fire or explosion endangering the user s of the device . One possible method of avoiding thermal runway is to change to composition of the battery by removing the liquid electrolyte material and rep lacing it with a solid - state fast lithium - ion conductor like that of the lithium garnet explored in this thesis. - ion batteries are made up of three major components , the electrodes o ne negative the other positive and the electrolyte ( Figure 1 ). These constit uents together The assembly of many electrochemical cells make what is colloquially known as a battery. The properties of the battery Figure 1 : Diagram of a simplified electrochemical cell in the charged state. The ion (green) is transported through the electrolyte, the electron ( yellow) travels through an external circuit, where it recombines with the ion in the cathode. 3 are thus dependent upon the specific reduction and oxidation reactions that occur bet ween the components within the cells. M uch of the work in the development of new battery technology is spent focusing on how to improve the performance of a battery by modifying the composition of the three major components in order to find an optimal comp osition. Operation of a battery requires the movement of both electrons and ion s from an electrode of higher potential to one of lower potential by means of an electrochemical gradient. The diffusion of ions within the gradient can be described as arising from Onsager like thermodynamics . 16 T he diffusion of mass and charge is coupled through the electronic potential of elections / holes and the chemical potential of a species or their vacancies. The potential is acted upon at the surface of an electrode where a redox half reaction occurs generatin g lithium - ion s and election s . T circuit . S imultaneously , an electron travels through a wire from one electrode to the other . At the opposing surface , an electron will recombine with the lithium - ion by a reduction half reaction completing the coupled mass - charge transfer. Importantly, the system as described implies a passive electrolyte with complete reversibility of chemical and physical reactions in the electrodes . That is to say, d uring cell cycling the electrolyte only transports ions between electrodes and does not directly partake in chemical reactions . The major challenge in designing an electrolyte is that a ny change in composition , volume , or texture at the electrode s must be completely accommodated by the electrolyte. Anyone who owns a cell phone can understand why this simple model is not quite true. Over the lifecycle of the battery it is quite common to experience a diminish ed c apacity , or loss of energy accessible in the battery . The major culprit being parasitic side reactions that convert electrolyte material and degrade electrodes as energy is transferred in the cell . 4 Having knowledge of the materials used currently in lithi um - ion batteries is therefore an important first step in determining what can be tweaked and modified to achieve better performance. Positive electrodes are typically crystalline solids dispersed in a film, taking the forms of Li x MO y or Li x M(PO 4 ) y , where M is a transition metal such as Co, Ni, Mn, or Fe. 17 Cobalt is most commonly used due to its superior performance and relative stability. However, ost to mine and refine result in a price premium. Using binary or ternary substitution at the M site can help compensate for the cost of cobalt by using cheaper m etals and in some cases e ven improve the specific capacity . Such is the case for lithium cobal t manganese oxide and lithium nickel cobalt aluminum oxide . In any case, typically Li/Li + half reaction potential is around 4 volts and specific capacity will range from 150 to 200 mAh g - 1 . On the other electrode, c ommercially available negative electrode s are commonly made from group IV materials: carbon, silicon, tin, germanium, or the crystalline solid lithium titanium oxide (LTO). 17 These can further be broken up into two main types, intercalation (C, LTO) and conversion materials (Si, Sn, Ge). Intercalation electrodes exhibit small lattice strain upon in sertion because of voids within the structure hav e enough space to accept a lithium . The most commonly used negative electrode material is graphit ic carbon, which allows lithium ion s to insert themselves between the sheets of layered carbon. combi nation of relative ly low cost, high quantity, high electrical / ionic conductivity/diffusivity, and low Li/Li + reaction potential make it a good choice as a negative electrode in most applications . The major drawbacks of carbon relate to its reactivity with some electrolyte materials which can result in exfoliation of graphite sheets leading to diminished capacity over time. While conversion electrodes exhibit much higher capacity, they operate in a manner in which the electrode undergoes chemical reactions. The major disadvantages being imperfect reversibility of the 5 reactions and extreme changes in electrode volume upon lithium o xidation/reduction which can diminish interfacial contact and isolate particles m aking them inert. The final component of an electrochemical cell to be discussed is t he electrolyte . I n most modern lithium - ion batteries , a nonaqueous solution of organic esters and carbonate solvents are used to dissolve a lithium salt. The salt will dissociate in solution releasing lithium - ion s which can be used to facilitat e the redox reaction s at the electrodes . A typical solution would be that of ethylene carbonate (EC) and dimethyl carbonate (DMC) with a salt such as LiPF 6 . 18 Here t he cyclic EC is used to dissolve the salt, while the linear DMC is used to lower the melting point of EC optimiz ing the transport of ions. It is common during the first cycling of the cell to form what is known as the solid electrolyte interface (SEI) 13, 19 . This complex material is a product of the electrolyte and electrode reacting to form a passivating layer. This passivatin g layer helps control the rate of electrode degradation but impedes lithium activity resulting in slower kinetics and loss of capacity. The exact formulation of the SEI is an active area of research and dependent upon the composition of both the electrolyt e and electrode partaking in the reaction. 19 - 20 The organic nature and liquid phase of modern electroly tes brings with it a plethora of problems. During normal operation, patristic side reactions between the electrolyte , SEI, and electrodes can result in capacity loss during cycling. Exceeding a safe operating voltage can cause the SEI to break down, fragme nt, and lead to further degradation of electrodes resulting in capacity loss. 13 The narrow range of thermal stability of the solvents requires supplemental cooling systems that add complexity and mass. 7 Manufacturing concerns like moisture contamination can result in hydrolysis reactions that break down the SEI. 8 Most critically, if any 6 one of the above systems failures arise, thermal runaway can occur igniting the organic solvent causing potential for injury of the end user. 21 solid - state batteries using a solid electrolyte are ag ain gaining traction in the literature. 22 - 23 In principal, the operation of the battery would remain the same, excep t the liquid electrolyte would be replaced by a solid fast - ion conductor, removing the hazardous organics. A wide range of solid electrolytes have been investigated exhibiting conductivities comparable to those of liquid systems above 10 - 3 S cm - 1 . 2 Figure 2 plots the conductivity for a range of solid electrolytes at various temperatures. The best performing of these, Li 10 GeP 2 S 12 (LGPS) , shows extremely promising room temper ature conductivity of around 10 - 2 S cm - 1 but suffers from moisture sensitivity and a stability window possibly as small as 1.7 to 2.1 V 24 . Similar concerns of Garnet LLZT Figure 2 : Conductivity of various electrolyte materials compared to that of LLZT 2 . The blue box designates the room temperature conductivity for the various materials 7 stability and shared with other high performing materials like the polymer/gel and phosphate - based electrolytes making oxide - based compounds an interesting route of study. Of the solid electrolytes in Figure 2 , the lithium garnet series Li 7 - x La 3 Zr 2 - x Ta x O 12 (LL ZT ) has been of great interest because of its high ionic conductivity, high thermal stabi lity, passivity toward common electrodes, and flexibility to incorporate various structural ion substitutions that allow for lithium concentration tuning. 25 - 33 Figure 2 places the conductivity of the garnet series about an order of magnitude lower than conventual electrolytes or LGPS at room temperature . Bridging the conductivity g ap is of interest to researchers and engineers because it will allow for the incorporation of solid electrolytes without diminished performance as compared to current batteries. The means for increasing the conductivity of LLZT has been focused on doping o r su bstituting aliovalent ions modifying the composition of the framework atoms and the distribution or concentration of lithium in the crystal. Experimental measurements have allowed scientist s to identify the structure and conductivity of the garnets, b ut to obtain a more comprehensive understanding of the mechanism of diffusion requires an atom ic - level precision . M aterial simulation is therefore a required technique to better understand how the movement of individual atoms gives rise to the macroscopic properties. Through simulation , we are capable of tracking both the position and velocity of every atom which allow s us to apply analysis techniques derived from statistical mechanics and observe their relation to properties at the thermodynamic limit . By simulating this material, we hope to obtain properties in a greement with experiment al values with respect to diffusivity, conductivity, and phase change phenomena . Once confirmed, we then take a deeper look at the distribution of species to better understand the role of the local environment of each species and how it modifies the dynamic properties . With a better understanding of how the 8 garnet system works , not only will we be able to provide suggestions o n how to improve the performance of this material, but also take our understanding of simulating materials and expand to other ionic conductors of interest. 1.2 The l ithium g arnet c rystal T he name garnet for these electrolyte materials originates from the mineral, Ca 3 Al 2 (SiO 4 ) 3 , which share s the same crystallographic space group s Ia3d and I4 1 /acd . The formulation A 3 B 2 C 3 O 12 is commonly used in geology and crystallography, representing the eight, six, and four - fold coordination with oxygen respectivly . 29 It is common among researchers studying the material to rearrange this formula to signify the importance of lithium to the form Li x A 3 B 2 O 12 , w here the lithium concentration is derived by the oxidation state of the A and B cations. W ith regards to the work presented the A cation is lanthanum in its 3 + state while B is either tantalum or zirconium in the 5 + and 4 + states respectively. The resulting formula of Li 7 - x La 3 Zr 2 - x Ta x O 12 (LLZT) can thus describe the range of compositions explored in this work. Figure 3 : Space filling model of L5LT (left), La (green) and Ta (brown) make a framework structure for lithium (blue) to diffuse. (right) Lithium tetrahedral (pink) and octahedral (blue) sites in the lithium garnet (Wang 35 ). (a) (b ) 9 G arnet was first discovered to be a potential solid electrolyte in 2003 with the work reported by Thangadurai, Kaack, and Weppner in 2003. 25 Ba sed upon the structural refinement of Li 5 La 3 M 2 O 12 crystals , where M is a metal of 5 + valency, it was s hown that chains of La 3 M 2 O 12 are surrounded by interconnected tetrahedral and octahedral sites with random lithium occupation. It is these dispersed vacancies that give rise to the potential of lithium diffusion network which is show in Figure 3 . The blue spheres in Figure 3 (a) show all the possible lithium positions within the crystal on the left, while the right figure shows the local spa t ial arrangement of the tetrahedral (pink) and octahedral (blue) sites. The lithium sites are commonly referred to by their Wyckoff positions in the cubic phase, the tetrahedral sites being the 24 d site, and the octahedral being the 48 g site. After the initial reporting of high ionic conductivity in the garnet class of materials , there was great interest in finding out what other cations could be substituted in the crystal in hopes of improving the conductivity an d reducing material cost . A review by Thangadurai et al 34 gives a thorough reporting of the various experimental results for approximately one hundred different garnet compositions . It is evident from the number of accommodating substitutions that the garnet structure provides a robust platform for controlling lithium concentration through multiple six coordination with oxygen forming octahedrons. Substituting tellurium with its 6 + oxidation state allows for a minimal concentration of lithium of three per formula unit (p.f.u). While lithium concentrations exceeding seven p.f.u. have been achieved by co - doping with the crystal with 3 + ions such yttrium or scandium with 4 + ions of zir conium, tin, or hafnium. 10 A maximal conductivity is reportedly achieved with lithium concentrations between six and seven 35 with dramatically lower room temperature conductivities for compositions at the end of the series w hen the [Li] = 3 & 7 34 . The dramatically lower conductivity arises from lithium ordering on the tetrahedral and octahedral sites that form the lithium sublattice . As shown in Figure 4 (a ) when [Li] = 3, lithium fully occupies the tetrahedral sites while leaving all octahedral sites unoccupied. It was shown through simulation work by Xu et al 36 that a large energy barrier exists within the octahedral that join tetrahedral sites making the conduction of lithium unfavorable. As the concentration of lithium is increased Figure 4 (b) , the only sites available for occupation are the octahedral sites which are seemingly unfavorable when compared to tetrahedral sites . The disorder in the lithium sublattice is hypothesized to arise from Li Li repulsion which results in a balancing of site preference to minimizing lithium interactions . 37 In this work an examination of how these forces balance is explored by applying a regular solution model to the disordered compositions, but much work is still to be done in understanding this compli cated relationship. Figure 4 : Diagrams of lithium occupancy for three different compositions of LLZT at 300 K . (a) Li 3 La 3 Te 2 O 12 , (b) Li 5 La 3 Ta 2 O 12 , and (c) Li 7 La 3 Zr 2 O 12 . (a) (b) (c) 11 For the [Li] = 7 compositions ordering again plays an important role in describing why the conductivity is so low compared to other compositions ( Figure 4 (c)) . In the case of [Li] = 7 structures, the crystal symmetry changes from Ia3d to I4 1 /acd at room temperature , resulting in a tetragonal crystal . The transformation results in a s plitting of the 24 d tetrahedral into two sites . A n 8 a site that has full lithium occupancy, and the 16 e site that has zero lithium occupancy. T he 48 g octahedral site also splits into 32 g and 16 f sites both with full occupancy. The defining characteristic separating the two octahedral Wyckoff positions being that the 32 g octahedra bridge 8 a and 16 e tetrahedra, while 16 f bridge two 16 e tetrahedra. A more general refinement splits each 48 g site into two 96 h sites which can help account for the displacement of lithium off the crystallographic ideal. 37 A simplified schematic of these symmetry positions as they relate to lithium are shown below in Figure 5 (a) and (b) with the net result of ordering in Figure 4 . Previous work by Wang et al. 37 provided structural refi nement data and analyzed the occupancy of lithium on its crystallographic sites using a combination of neutron diffraction and energy minimization simulations. Essentially, they showed that there appears to be a preference for tetrahedral occupancy at low lithium concentrations an d temperatures due to the increased symmetry of the site . Lithium o ccupancy into octahedral sites then results as i s a means to Figure 5 : Schematic showing the relationship between possible lithium Wyc koff positions in the cubic (left) and tetragonal (right) crystals . Filled circles are occupied lithium sites and hollow circles are unoccupied sites. 24 d 48 g 24 d 96 h 96 h 8 a 16 e 16 f 32 g (a) (b ) 12 minimize the Li Li interactions by populating the larger and less coordinated octahedral sites. This work provide d a good understanding of how one can expect lithium to act under static conditions and is to be expanded upon in the proceeding work of this thesis as to better understand the dynamic properties of the garnet system . Properties such as how the preference f or certain Li Li configurations change with composition and temperature, how these connecting the tetrahedral and octahedral faces, and how changes in the lith ium sublattice relate to the changes in the overall structure and phase of the garnet crystal. 1.3 Tetragonal to cubic transformation in Li 7 La 3 Zr 2 O 12 1.3.1 Experimental observations T he f irst reporting of Li 7 TLa 3 Zr 2 O 12 ( L7LZ ) compositions showed a cubic phase with a room temperature conductivity of 0.5 mS /cm . 38 It was later discovered that Al 3+ ions leeched from the crucibles during the synthesis process , replacing tetrahedral lithium , resulting in the cubic phase the high conductivity . 30, 39 - 41 Subsequent experimental work reported tetragonal phase garnet with orders of magnitude lower conductivity. Groups followed up on this work by purposefully doping Al 3+ or Ga 3+ ions which causes disorder on the lithium sublattice by repelling neighboring octahedral lithium and removing three total lithium in order to maintain charge neutrality. Electron diffraction studies performed by Thompson et al. 42 establish that a critical concentration of [Li]=6.5 is necessary to achieve a room temperature cubic phase with high ionic conductivity. It is proposed that a similar critical concent ration is needed to achieve a room temperature cubic phase trough tantalum doping on the zirconium site. W hile experimental results suggest that this is indeed the case, an in - depth investigation into the mechanisms controlling this transformation through modeling warrant exploration. 13 P ure L7LZ upon sufficient heating will transform from the tetragonal phase to a cubic one exhibiting high conductivity . The phase transformation is thought to be a first order phase transformation expressing volume conservati on , with the added energy contributing to an increase in entropy attributed to lithium disordering . 43 The cause of the transformation is not explicitly known and is a major point of interest for thi s thesis work. Matsui et al 4 report s the phase trans formation temperature to be between 600 and 650 o C through high - temperature X - ray diffraction ( XRD ) and impedance spectroscopy. Larraz et al. 43 report a temperature around 650 C in their own XRD and differential scanning calorimetry (DSC) experiment. While Wang et al. 1 identify a transformation around 610 - 630 C depending upon the he ating or cooling of the system via impedance spectroscopy . These reported temperatures are fairly consistent with pure L7LZ as will be shown in the simulation work presented in the following sections. 1. 3.2 Simulations of garnet in literature S imulating th e transformation from tetragonal to cubi c is therefore an important reference point for model verification. Bernstein et al 44 approached modeling the phase transformation using molecular dynamics derived from density functional theory (DFT) calculations. They determined vacancy free L7LZ would transform between 800 K and 1000 K. Introducing lithium vacancies while keeping zirconium concentration constant resulted in a depressing of the transition temperature to below 300 K when the critical concentration was o btained. They postulate that the transformation to a tetragonal symmetry occurs in order to preserve ZrO 6 octahedra orientation due to increased columbic interactions between Li Li pairs in the ordered state. While this methodology fits the narrative that lithium vacancies are the reason for depressing the tetragonal to cubic transformation, their model does not account for cooperative 14 effects tantalum may play in the phase transformation witch its smaller ionic radius and higher formal charge . Other MD - D FT studies on garnet materials have been performed , namely the work by Xu et al 36 , Miara et al 5 , Meier et al 45 , Jalem et al 46 - 48 , and Rettenwander et al . 49 In general these studies avoid making thei r models transform from one phase to the other by simulating temperatures above the transition point or by explicitly making two crystals, one for the ordered tetragonal phase and one for the cubic phase. This is done to focus more clearly on the location of dopants and the motion of lithium in the crystals versus the transformational effects as explored by Bernstein . Each of these studies only employ one type of electron exchange - correlation (XC) energy functional, leaving to question how optimal each simulation was with respect to modeling the complex bonding environment in these crystals. This question is addressed in Chapter 4 of this thesis in my work looking solely at the role of the XC functional in modeling the L7LZ high temperature phase transformation. Minimum energy simulations for the Ta/Zr series of garnet performed by Wang 37 et al using General Utility Lattice Program (GULP) showed how the configurational energy of static structures of garnet distribute with varying levels of lithium concentration. In this experiment 5000 initial structures wit h random Li, Ta, and Zr distribution were optimized for each composition. T he cubic garnets with [Li] = 5 - 6.25, exhibit mostly gaussian distributions and inspection of their atomic positions suggests all stable disordered structures. As the concentration o f lithium is increased further to 6.5 - 7, the single gaussian begins to split into two peaks, with the lower energy peak becoming increasingly spiky as more lithium is added. Structures that optimized in this spike all exhibit identical structures upon opti mization which are associated with the ordered tetragonal phase. Importantly to the work presented in this thesis, increasing the 15 size of simulation from a 1 x 1 x 1 unit cell to a 2 x 2 x 2 cell gave statistics consistent with classical statistical thermo dynamics and provided a solid foundation for evaluating the dynamic properties of these materials using molecular dynamics. The first attempt at modeling L7LZ using classical molecular dynamics was done by Adam and Rao 50 whose model was fit to diffraction data that supported a much lower phase transition temperature of only 400 K . The model is capable of identifying possible diffusion pathways, but values of activation energy, conductivity, and site occupancy may not be reliable due to the lower than predicted phase transition . As will be seen in the proceeding work of this thesis , getting a single model to accurately capture the properties of these materials is complex. Initial work by us in Wang et al 51 and Klenk et al 52 modeling L5LT and L7LZ respectively, using classical molecular dynamics with a core - shell model, are able to reproduce lithium dynamics in agreement with experimental quasi - elastic neutron scattering experiments for L5LT . 53 But, when trying to predict accurate lattice volumes fo r L7LZ , the model over binds the atoms predicting a smaller volume and distorted lattice. The following work in Klenk et al 16 using a core only MD model, is much better at approximating the transformation from tetragonal to cubic in L7LZ , and present s compelling structural information for intermediary compositions, though the dynamics of lithium are likely overestimated. This is of course the difficulty in modeling real materials , it is important to identify what they perform well and if it is possible to generate useful information with them. From the work highlighted above you can see that simu lations of the garnet can be divided into two main categories. First there are the DFT studies that use the higher - level calculations to investigate individual atomic interactions in order to better understand the energetics of diffusion pathways and the e ffects of dopants on the arrangement of atoms . The 16 other set being classical molecular dynamic studies that take advantage of the more efficient calculations in order to run larger and longer simulations . This allows for better statistical sampling and ultimately hopes to better approximate macroscopic properties . 1.4 Experimental measurement of c onductivity and d iffusivity 1.4.1 Conductivity I n order to measure the conductivity of a solid electrolyte like the lithium garnet one would most commonly use electro chemical impedance spectroscopy (EIS) . Typically, an electrical stimulus is applied to a system , either a known voltage or current , and b ased upon the response of the system , it is possible to deduce what kinds of electrical processes are occur. When a sinusoidal input is applied, like AC voltage, to a linear system, the output current will have an identical angular frequency. Equation 1 The impedance Z is comprised of two components one real (Z R ) and one imaginary (Z I ) . Equation 2 Plotting the real vs imaginary parts of Z( ) in the complex plane will give a characteristic semi - circle when applied to an ionic conductor such as the lithium garnet. By defining an equivalent circuit, typically a constate phase element /resistor in parallel with a capacitor , it is possible to fit the measured spectrum to the model circuit and ide ntify the bulk, grain boundary, or total resistance of the system. Taking the geometry of the system into account and obta ining the dimensionless resistivity , then taking the inverse gives the conductivity of the sample. Because 17 collecting the conductivity is a fairly easy experiment to perform there are many reports of the conductivity for a verity of compositions. An exhaustive list of measured conductivities can be found in the review by Than gadurai. 34 It is assumed that the garnet series follows Arrhenius behavior as the temperature increases, and fitting the conductivity vs inverse tempera ture is used to obtain activation energy values for the solid electrolyte. 1.4. 2 Di ffusivity U nderstanding the Li diffusion mechanism in this material is essential to optimizing Li transport behavior and realizing its full potential as a solid electrolyte. Experimentally, it is very difficult to directly measure the diffusivity and requires the us e of exotic techniques such as nuclear magnetic resonance (NMR) 31, 54 - 58 and muon - spin relaxation ( µSR ) 59 - 60 C onventional powder diffraction analysis 3 has been combined with maximum entropy methods to study the Li distribution in lithium garnet oxides 61 . This work found a channel of lithium that delocalizes along 24 d 96 h 48 g 96 h 24 d site chains, through both tetrahedral and octahedral cages. Although giving insight into the lithium diffusion pathway, this experimental approach does not directly measure lithium dynamics. The translational dynamics of Li in lithium garnet oxides have been the subject of intense investigations using nuclear magnetic resonance (NMR) 31, 54 - 58 and muon - spin relaxation ( µSR ) 59 - 60 measurements. For example, the residence time of Li at crystallographic sites was extracted by applying the Bloembergen - Purcell - Pound (BPP) equation to the longitudinal/spin - lattice relaxation (SLR) time 31, 54 - 55 . Alternatively, Li self - diffusivity was obtained by applying the Stejskal - Tanner equation to pulse gradient spin echo (PGSE) NMR measurements 56 - 58 . Similarl y, only the residence time of Li could by extracted from the dynamic Gaussian Kubo Toyabe function by µSR measurements 59 - 60 . If we consider that diffusivity ( D ) is related to both 18 length ( l ) and residence time ( t ), such as by , SLR - NMR and µSR techniques are extremely limited in that they only yield information for , and PGSE - NMR only provides information at the long - range limit of . To probe diffusion both spatially and temporally, scattering - based techniques such as qu asi - elastic neutron scattering (QENS) are essential. QENS measures neutron scattering intensity as a function of momentum ( Q ) and energy ( E ) transfer (or Fourier time in a spin - echo measurement) of neutrons interacting directly with atomic nuclei. In terms of atomic dynamics, its result s contain the neutron - weighted self - part of the van Hove correlation function 62 , i.e. the probability of findin g a particle after it travels a length of r for a time duration of t , shedding light on self - diffusion. Before this work was initiated there has been only one QENS report on lithium garnet oxides 59 in which only the elastic intensity of the measured spectrum was quantitatively treated and the spatial and temporal features were not discussed. The goal of this thesis is to tie together simulation and experimental approaches in order to present a compelling narrative for how the dynamic properties of the lithium garnet series are related to the crystalline structure with respect to lithium concentration and arrangement . Further, a n investigation into the effects of Ta/Zr dopin g using simulations over a wide range of concentrations ha s not been reported and currently there is no verified explanation as to why the we see a maximum conductivity at the concentrations we do . The procedures and analysis techniques utilized in this work are of great interest to those specializing in simulating materials an d are generalizable to other materials. 19 CHAPTER 2 : Materials and m ethods 2 .1 Sample preparation L i 5 La 3 Ta 2 O 12 p owder was synth esized by conventional solid - state reaction methods. 7 LiOH · H 2 O (Cambridge Isotope Laboratories) and La 2 O 3 (Sigma - Aldrich) were dried at 215 o C and 600 o C, respectively before mixing. 7 LiOH · H 2 O, La 2 O 3 , ZrO 2 (Alfa - Aesar), and Ta 2 O 5 (Alfa - Aesar ) were added in stoichiometric amounts with a 5% excess of 7 LiOH · H 2 O to account for lithium loss during synthesis. Sample precursors were mixed in polyethylene jars and ball - milled for 24 h in isopropanol with yttria - stabilized - zirconia grinding balls. The dried mixture was placed in MgO crucibles and fired at 900 o C for 10 h. Upon completion of the synthesis , powders were placed in an argon atmosphere and vacuum packaged in metallized bags for storage to minimize moisture contamination. The powders were ex amined using X - ray diffraction with a Bruker D8 Advance diffractometer and shown to be phase - pure. 2 .2 Quasi - elastic neutron scattering Q ENS data were collected on the backscattering spectrometer (BASIS) 63 at the spallation neutron source (SNS) at Oak Ridge National Laboratory (ORNL). Approximately 10 g of powder was sealed in an annular Al container and was initially heated to 700 K to remove any moisture and co ntaminating gas accumulated during handling. Spectra were collected in 100 K decrements to 300 K and then at 50 K to define the resolution function. Data were converted from time - of - flight to energy transfer spectra at selected Q using Mantid software package 64 . S( Q , E 100 to Q over 0.1 Å - 1 from 0.25 Å - 1 to 1.95 Å - 1 . The QENS analysis was performed using the Data Analysis and Visualization Environment (DAVE) package 65 . 20 2 .3 Molecular dynamic simulations T he majority of the work presented in this dissertation is the result of molecular dynamic (MD) simulations comprising of three different simulation models for evaluation. Two of which are pair - potential models, and one that uses DFT. Classical MD, being a simpler technique is based upon the idea that only the position of the particles is necessary to calculate the potential energy of a system, while DFT requires solving for the potential energy of a system based upon the electron density. that , with F being a force, m, the mass, r the position, t the time, and U the potential energy. From this, if you can predict the next position of an atom after some small time differential it should be possible to calculate the force experienced by each particle. To calculate the potential energy (U) of a system as a summation of pair - wise interactions where . In this work we u se a pair - potential model that simulates the long - range Coulombic interactions and short - range Buckingham interactions as Equation 3 where A, , and C are fitting parameters, r ij is the distance between species i and j, z is the charge, and 0 is vacuum permittivity. Once the force on each atom is calculated dividing by its mass will give the acceleration experienced by each atom setting up the possibility solve for the atoms next position through Verlet integration. Equation 4 21 This process of calculating the potential energy, to then calculate the next velocity and position is repeated until a specified number of steps has been performed. 2 .3.1 Core - shell model T he first model to be discussed is the Dick - - polarizable species like Ta and O. A simple harmonic oscillator is then used to connect the cores describe the bonding environment within the garnet at the expense of shortening the time - step used in the integration. Our Buckingham, electrostatic, and core - shell parameters for the core - shell model can be found in Table 1 . Initial average atomic position and lattice parameters for the garnet series were calculated from neutron scattering experiments performed by Wang et at. 37 From these structures, we randomized Li, Ta, and Zr distributions and performed static energy minimization at 0 K using the General Utility Lattice Program (GULP). The majority of the model terms A, , and C, were taken from literature, while the preexponential terms (A) in Table 1 for Zr - O and Ta - O were obtained from potential parameter fitting from neutron diffraction structures at 10 K. The energy minimization of structures was performed in two steps, first to relax only the position of the electron shells, and then a relaxation of all species, to improve convergence to the local minimum. The minimization utilized conjugant gradient then hessian inversion for solving for the next local minimum when the difference in the gradient was above or below 2 eV respective ly. 22 It is important to note that when dealing with compositions of LLZT with a [Li] < 7 it is computationally infeasible to obtain an absolute minimum energy due to the disordered nature of the lithium sublattice. There are between 40 and 56 lithium atoms per unit cell to be distributed about 72 possible ideal sites meaning the number of permutations of lithium occupancy is immense. In order to obtain a representative structure of a unit cell in the bulk, 1000 randomized 2 x 2 x 2 simulation cells were gene rated and optimized. The energy of the optimized structures was placed in a histogram and a structure was taken from the most probable energy. MD simulations were performed with the DL_POLY Classic software package 66 at temperatures ranging from 300 K to 1400 K . To calculate the trajectories of the ions, the vel ocity V erlet integration method was employed with a 0.25 fs time - step. The first step in our simulation process was to take one of the optimized structures from the GULP evaluation and heat it to 1400 K for L7LZ, and 1200 K for L5LT, under constant number, stress/pressure, and Table 1 : Core - shell Buckingham and spring pair interaction parameters for LLZT 23 temperature (N T/NPT). The Berendsen thermostat and barostat was used with relaxation time constants of 0.1 and 0.5 ps respectively and a stress/pressure of 0 GPa. We then sequentially cooled the resulting structure to 300 K in 100 de gree decrements for L7LZ, and at temperatures of 1100, 900, 700, 550, 475, 350, and 300 K for L5LT. The average lattice parameters were taken from the last 50 ps of the NPT simulation. The final configurations of the NPT simulation were then rescaled to th e calculated average lattice parameters for use in constant number, volume, and energy (NVE) simulations. NVE simulation times varied depending upon the particular study, in Wang et al. 51 (L5LT) and Klenk et al. 52 (L7LZ) 100 ps of equilibration and 1ns of pro duction was performed with trajectory output every 0.1 ps. An additional 2.5 ps of production simulation with output every 2.5 fs was performed at 300, 1100 K (L7LZ) and 550 K (L5LT) to better refine more subtle details of diffusion. In Klenk et al. 53 the total number of temperatures inve stigated was reduced to 400, 475, 550, 700, 900, and 1000 K. The same equilibration was used as above but the production run time was reduced to 500 ps for 550, 700, 900, and 1000 K. For 475 K a production time of 5 ns was used while at 400 K the run was e xtended to 10 ns. The mass of the O shell was increased to 0.4 AU and the time step for simulation was reduced to 0.1 fs at 400 K in order to keep the system adiabatic. 2 .3.2 Core - only m odel C omputational constraints with DL_POLY Classic restricted our simulations to a maximum size of 2 x 2 x 2 due to the parallelization using the replicated data strategy. Here the positions and velocities are stored on every processor increasing causing exce ss memory usage. We saw no appreciable improvement in running DL_POLY in parallel and choose to limit our 24 simulations to run on a single processor. A more efficient means of allocating memory in a simulation would be to employ a domain decomposition strate gy that separate the storage and computation of forces into subdomains determined by the spatial representation of the unit cell, and the number of processors used for simulation . We choose to use the Large - scale Atomic/Molecular Massively Parallel Simulat or ( LAMMPS ) 67 as a means of increasing our simulation size in order to explore finite sizing effects, and the effects of Ta/Zr substitution for compositions with [Li]= 5, 5.5, 6, 6.25, 6.5, 6.75, &7. We were able to perform simulations with up to a size of 4 x 4 x 4 for L7LZ, and as will be shown 3x 3 x 3 simulations cells were deemed large enough to observe convergent behavior for all other composition. Due to imitations in LAMMPS not allowing for the use of our p revious core - shell model we instead employed a partial charge model that is based on first - principles density functional theory (DFT) calculations. Specifically, we obtained the partial charges of atoms in the DDEC (Density Derived Electrostatic and Chemi cal charge) scheme 68 , from results of DFT calculation with the VASP package 69 , 70 , 71 , 72 em ploying the Projector Augmented - Wave (PAW) method. 73 , 74 The simulation cell was a 1×1×1 supercell of Li 7 La 3 Zr 2 O 12 with 192 atoms. Structural parameters were taken from experimental values in literature. 37 A self - consistent calculation was performed with the Generalized gradient approximation (GGA) of Perdew - Burke - Ernzerhof parametrization (PBE) 75 , 76 as the exchange - correlation functi onal, a 600 eV cut - off energy, and at the Gamma point. Valence electron configurations for La, Li, O, and Zr atoms are 5s 2 5p 6 5d 1 6s 2 , 1s 2 2s 1 , 2s 2 2p 4 and 4s 2 4p 6 4d 2 5s 2 , respectively. The DDEC charges of all atoms in the same species were averaged to yield v alues of 2.53 (La), 1.01 (Li), - 1.70 (O), and 2.86 (Zr). In order to facilitate the investigation of Li 7 La 3 Zr 2 O 12 - Li 5 La 3 Ta 2 O 12 we slightly offset these DDEC charges to values reported in Table 2 . The interaction parameters were refined against 25 the experimental structure 37 by th e GULP package 77 , 78 , with initial values taken from the reference. 48 The force - field parameters are listed in Table 2 Table 2 : Core - only force field parameters Buckingham Pair Interaction with Oxygen Species Mass z ( e ) A ( e V) (Å) C ( e V Å 6 ) Li 6.9400 1.00 1087.29 0.260 0.00 La 138.9050 2.50 2075.26 0.326 23.25 Ta 180.9479 3.65 1177.30 0.322 0.532 Zr 91.2240 2.65 1650.32 0.311 5.10 O 15.9990 - 1.65 4870.00 0.267 77.00 For the disordered compositions of LLZT, 200 3 x 3 x 3 simulation cells were generated with randomly dispersed Li, Ta, and Zr. The local energy minimization was obtained using GULP static energy minimization based upon structures from neutron diffraction. 37 Due to the increased size of the simulations, it was no longer efficient to use Hessian inversion in order to obtain the gradient in energy for these structures. Instead, only the conjugant gradient method was used to offset the atoms from their ideal pos itions to obtain convergence within 0.001 eV per structure. 2 .4 Finite s ize effects in L7LZ I nitial tetragonal L7LZ simulation cells consisting of 1×1×1 (192 atoms), 2×2×2 (1536 atoms), 3×3×3 (5184 atoms), and 4×4×4 unit cells (12288 atoms) were generated using the same structure as the DFT calculation. The Ewald summation of long - range interactions was calculated using a cutoff radius of 6 Å for 1×1×1 cells and 12 Å for 2×2×2 and larger cells. Anisotropic NPT simulations were carried out by heating the i nitial tetragonal structure to 1400 K and sequentially cooling the simulation cells in 100 K decrements until 300 K was achieved. 26 For the 3×3×3 and 4×4×4 cells, we implemented additional sequential heating and cooling from 800 - 1000 K in 25 K increments to study the hysteresis. NPT simulations ran for 75 ps with a 1.0 fs time - step employing the Nose - Hoover thermostat and barostat 79 - 82 with relaxation time constants of 0.05 p s and 0.25 ps respectively. Lattice parameters were extracted from the last 50 ps. Once average lattice parameters were determined, NVE simulations were performed for 1 ns utilizing the same 1.0 fs time - step. The atomic trajectories are saved every 10 f s. 2 .5 Effect of Ta/Zr substitution I n accordance to the procedure used for the core - shell model outline above, 3 x 3 x 3 LLZT compositions with [Li]= 5, 5.5, 6, 6.25, 6.5, & 6.75 were generated by selecting a structure from the most probable energy in a histogram of energy distribution. These structures were simulated under anisotropic NPT isotropic conditions at 1400 K and sequentially cooled in 100 degree decrements down to 300 K. As with L7LZ these simulations ran for 75 ps with the same time - step, bar ostat, and thermostat settings. The average lattice parameters were taken from the last 50 ps from anisotropic NPT calculations, and averaged to form cubic lattices. The atomic positions from the last frame of the NPT simulation was rescaled to the average lattice parameters calculated to be used in NVE simulations. Initial velocities were generated randomly based upon a gaussian distribution, and equilibration was run for 100 ps. A time - step of 1 fs was used to calculated trajectories for up to 1 ns, with data output every 10 fs. 2 .6 Density f unctional t heory: XC f unctional analysis 2 .6.1 Background R esearch focusing on simulating physical and chemical proper - ties of materials has grown tremendously as the raw computational power and algorithm efficiency have improved drastically over the past decades. Density functional theory (DFT) 83 - 84 has resultantly become a 27 common simu lation technique to predict the thermodynamic and kinetic properties of a variety of simple and complex materials based upon the electron density (r) . Starting with the Hamiltonian ( al energy (E) under the Born - Oppenheimer approximation is Equation 5 where is the kinetic energy of electrons, is the electron - electron interaction, is the electron - nuclei interaction, and is the interaction potential between nuclei. Kohn - Sham takes this calculation of energy and puts in in terms of the electron density Equation 6 with the kinetic energy external potential , and Hartree energy known, the electron exchange - correlation energy (XC hereafter) is an unk nown parameter. The XC functional is comprised of two parts one pertaining to the electron exchange energy (F x (s)), or the energy caused by the spatial separation of two elections with the same spin, and the correlation energy (F c (s)), the interaction energy of electrons from Coulomb repulsion. These parameters are presented as functions of s, the reduced density gradient. Equation 7 Within the DFT framework, the electron XC functional is not known a priori and approximate functionals have to be employed. This work takes functionals from the LibXC 85 package an d currently contains about 180 functionals with the expectation that more will added. 28 The simplest XC functional is the local density approximation (LDA) 84 which approximates the XC energy density at a local position by the value of a uniform electron gas 86 are generalized gradient approximation 75 (GGA) which includes the electron density gradient and meta - GGA which further includes the Laplacian of the density. Within each family of approximation, many different functi onal forms on the density, gradient, or Laplacian have been proposed based on the electron gas with slowly - varying density 75 , Airy gas 87 , or semiclassical neutral atom 88 . Due to the wide range of structure features and bonding environment within complex materials, choosing the proper functional for a system can be difficult. This is one of the reasons that much effort has been devo ted to test a variety of XC functionals for different solids to see how simulation results are close to experimental values (lattice constants, bulk modulus, cohesive energy, etc.), e.g. by Haas et al. 89 - 90 , Csonka et al. 91 , Hao et al. 92 , Labat et al. 93 , He et al. 94 , Rasander et al. 95 , Tran et al. 96 . However, most of these studies are focused on the static structure (0 K) of simple metals and binary compounds. Fast - ion conductors on the other hand are complex materials with at least one mobile species and of ten exhibit a phase transformation from a low - symmetry structure to a high - symmetry structure upon heating. At present, the effect of XC functionals on simulated materials properties of fast - ion conductors remains elusive. The objective of the present work is to evaluate a collection of XC functionals on the phase transformation of a model fast - ion conductor: lithium garnet oxide Li 7 La 3 Zr 2 O 12 ( L7LZ ). L7LZ was chosen because of its high lithium - ion ic conductivity, importance to solid - state battery research, and characteristic phase transformation from a tetragonal to a cubic phase around 900 K 38, 97 - 98 . L7LZ has been the subject to several DFT studies, e.g. Xu et al. 36 , Bernstein et al. 44 , Jalem et al. 47 , Meier et al. 45 , 29 Rettenwander et al. 49 , Miara et al. 5 , Kang et al. 99 b ut only one XC functional was employed in each of these investigations. In this work, we selected 14 different XC functionals (LDA, 13 GGA including 2 dispersion corrected ones) and recorded the dynamic response of the L7LZ system over the course of a simu lation to categorize if or what type of transformation occurs. We then compared the average crystal structure at different temperatures to experimental results in order to determine the best XC functional for this system. It is expected such investigation could also provide insight into studies of other fast - ion conductors. Table 3 : Exchange and correlation enhancement factors Name Exchange F x (s) Correlation F c (s) PBE 75 PBEc PBE2 100 PBEc PBEsol 101 PBEc with b=0.046 RPBE 102 PBEc SOGGA 103 PBEc RGE2 104 PBEc with b=0.053 WC06 105 PBEc 30 Table 3 (cont d) PBE - FE 106 PBEc with b=0.043 VMT84 107 PBEc PBE = 0.219516, GE = 10/81, FE = 0.346, PBE = 0.804, EL = 0.552, FE = 0.473 C WC06 = 0.00079325, VMT84 = 0.000023 2 .6.2 Simulation method T he first - principles DFT - based molecular dynamics computation was performed by the Quickstep code 108 implemented in the cp2k 109 package, which employs mixed Gaussian and pl ane wave (GPW) basis sets. The valence electron configurations were La (5s 2 5p 6 6s 2 5d 1 ), Li (1s 2 2s 1 ), Zr (4s 2 4p 6 5s 2 4d 2 ), and O (2s 2 2p 4 ) with Goedecker Teter Hutter (GTH) norm - conserving pseudopotentials. 110 - 111 The plane wave cutoff was 400 Ry and the Gaussian basis sets were molecular optimized double zeta - valence basis sets with a polarization function (DZVP) 112 . The simulation cell was a unit cell of Li 7 La 3 Zr 2 O 12 with 192 atoms and 1200 valence electrons, sampled at the gamma point. The constant number/pressure/temperature (NPT) simulation was performed for 6 or 12 ps (1 fs time step) at different temperatures and zero pressure, with a canonical sampling through v elocity rescaling (CSVR) thermostat 113 (time constant of 5 fs) and Martyna - Tuckerman - Tobias - Klein (MTTK) barostat 114 (time constant of 300 fs). Local density approximation (LDA) and thirteen generalized gradient approximation (GGA) XC functionals, 31 implemented in the original and modified LibXC 85 package were investigated. GGA functionals, except PW91 115 and AM05 87 , have 9 different PBE - like functional forms. The dependence of exchange enhancement factor F X (s) of these 9 functionals on the reduced density gradient, i.e. s, is presented in Table 3 They are mainly composed of two parameters, and . The parameter determines the behavior at small s, whose value was set to PBE in the PBE approximation by cancelling second - order exchange and correlation terms. Other values inspired by the gradient expansion of electron gas with slowly - varying density ( GE ) 101 or formation - energy fitting ( FE ) 89 whose value wa PBE to satisfy the Lieb Oxford bound 75 . Other tightened values EL ) 116 or formation - ene FE ) have been used 89 . The functional form of correlation enhancemen t factor F C (s) can be found in the literature 75 . Modified parameters of in F C (s) and other parameters in F X (s) are also presented in Table 3 . Finally, the dispersion correction in the Grimme D3(BJ) formalism 117 - 118 was applied to two of PBE - like functionals, PBE and PBEsol. 32 CHAPTER 3 : Molecular d ynamics a nalysis m ethods 3 .1 Diffusivity and c onductivity calculation methods T he efficiency of an electrolyte is measured through two main properties; diffusivity and conductivity. Ideally the two have an equivalency as described through the Nernst - Einstein relation in Equation 8 . Equation 8 Where is the conductivity, z i e charge of species i, D the diffusivity, is the nominal number density, k B is Boltzmann constant, and T the absolute temperature. However , Equation 8 is only true in the limiting case of dilute solutions of noninteracting particles though it is commonly applied to electrolyte materials like the garnet . A correction term donated in this work by - 1 , the inverse thermodynamic factor, can be applied to the right - hand side of Equation 8 to address the nonideal case as expanded upon later. This work takes the approa ch of separating the - diffusivity ( ; which captures the motion of individual particles, and the Fickian ( or collective diffusivity , that describes how the ensemble mass flux of particles moves . It is the ladder that can therefore be used make the mass - charge equivalency possible because it more closely relates to the drift velocity of the ions . To obtain the kinetic/dynamic properties for the lithium garnet series we use two main approaches . First , we can calculate the mean squared displacement (MSD) Equation 9 of ions at the long - time limit, where we expec ted to observe Fickian like diffusion as in Equation 9 . Where r(t) 0 is the initial position of a particle, r(t) t is the position at time t, n is the 33 dimensionality of the system , is the self - diffusivity , and t is the time . As seen in Figure 6 for L7LZ , the MSD is not perfectly linear and each line fitting will incur some error. The causes for the nonlinear behavior can stem from two caus es . E ither the length of time for simulation was not long enough , which results in too few sampling events for statistical evaluation , or the simulation cell is too small and we are observing errors incurred from the periodic boundary condition . T he nonlinearity was accounted for in this work by finding a region that is well into the simulation time and performing linear regression over a consistent time window. This approach however fails to address the range of possible diffusivities that could be r eported at each temperature. Taking 1200 K ( dark blue) as an example, if you were to fit the lines from 400 ps to 600 ps and 600 to 800 ps, the diffusivity calculated from the first set would be almost half that of the second set. Such lack of precision with this approach requires caution and is wholly dependent upon the proper selection of a time window. This problem is exacerbated when performing Figure 6 : Mean squared displacement of lithium ions in L7LZ at various temperatures during a 1000 ps simulation. 34 simulations at lower temperatures where the mobility of ions decays exponentially and jump times approach the length of the simulation. This simpler model of calculating diffusion was used predominantly in the early work by Wang 51 and Klenk 16, 51 - 52 and converted to conductivity via the m odified Nernst - Einstein equation ( Equation 8 ) . The limitations brought by obtaining the diffusivity via the MSD led to the development of a second approach, calculat ing the diffusivity and conductivity by way of the dynamic structure factor S(Q, ) as you approach the hydrodynamic limit , i.e. Q 0 and 0. 16 Practically for the sake of this work, we extrapolate from the smallest Q (Å - 1 ) for a given simulation size and equate that to be approximately the value at t he limit. A formal derivation of this method is outlined thusly. 3 .2 Calculation of collective diffusivity S tarting with 62 : Equation 10 where the flux ( J ) of species is the sum of the kinetic coefficient L multiplied by the gradient of the chemical potential of species . 119 With respect to the concentration gradient Equation 10 can be rewritten to take the form: Equation 11 with the term being the thermodynamic factor, T the absolute temperature, and k B the Boltzmann constant. 16 In the case of a charged species , the equilibrium chemical potential of Equation 10 can be expressed as where z is the charge of species , e is the 35 elementary charge, and is the electrical potential. Rewriting Equation 10 in these terms , the current flux can then be defined as: Equation 12 with being the species - species electrical conductivity of and . In the case of the garnet system there is only one mobile species lithium, so contributions from other atoms are generally ignored and reduces to Equation 13 . Equation 13 The kinetic coefficient L of Equation 10 , 11, 12 , & 13 is found by taking the species flux density based on the density autocorrelation at the hydrodynamic limit as Q goes to zero: Equation 14 Equation 15 where is the velocity and position of the k th atom of at time t, frequency, and V the volume of the simulation cell. The term represents the coherent scattering of species at time t and is the reformulation of the Kubo formula in Fourier space . 120 Equation 16 36 T he additional term w a is a weighting factor based upon whether we are looking at collective particle, charge, or mass correlation The correlation of velocity for species at time zero , with that of species at time t , (Q,t) , can take two forms as is expressed in Equation 17 , Equation 17 with being the longitudinal and transverse velocity correlations. Assuming isotropic diffusion behavior, a t the limit as Q 0 the longitudinal component goes more quickly to zero, while the relativity slower decaying transverse component can be used to obtain the kinetic coefficient by taking the Laplace transform in the time domain . 62 The rate of decay for the longitudinal and transverse velocity correlations is displayed for L7LZ at 1400 K in Figure 7 . The different weighting factors and their respective properties derived by the coherent correlation are shown in Table 4 . Table 4 : Properties derived by changing the weighting f actor in Equation 16 Property Weighting factor Equation Diffusivity (unity) Conductivity (charge) Sound/shear viscosity (mass) 37 3 .3 Practical considerations for fast - ion conductors W ith only one mobile species it is possible to rewrite Equation 11 law looking at the number flux for a collection of a single species. ( Equation 18 ) Equation 18 The terms between the summation and concentration gradient can then be considered as the Fickian or collective diffusivity . Using from Equation 18 it is possible to rewrite Equation 13 to the thermodynamically corrected Nernst - Einstein equation as in Equation 19 . Equation 19 Figure 7 : The Li Li velocity autocorrelation in the longitudinal (Black) and transverse (red) scattering directions. 38 Where is the conduction diffusivity, or the thermodynamically corrected collective diffusivity, which corrects for the finite scale of our small simulations cells and takes them to the thermodynamic limit of infinite media. This is in contrast to t he tracer or self - diffusivity ( ) which is commonly employed to make conductivity calcu lations . To obtain you can use MSD as described above or the velocity autocorrelation function for individual particles as in Equation 20 . Equation 20 Equation 20 is a generalization of Equation 9 defining the average displacement of all individual particles over some time. A simple analogy can be made to see the contrasting nature of finds how far individual pa rticles travel and then averages them to get an average displacement , while averages the position of the particles and measures how far they collectively displace or the rate at which the collection moves . We choose to replace in o ur calculations because the collective motion captured in is more representative of the current flux which is measured experimentally using impedance spectroscopy. We have outlined here three different ways of calculating the conductivity of lithium in the garnet system. The simplest and quickest route would be to calculate the self - diffusivity via the mean squared displacement, but the distinction between self and collective diffusivity is lost potentially making the analysis of conductivity unrepresentative of the true system dynamics . Rather, we could calculate the collective diffusivity by applying a weighting factor of one to Equation 16 , and proceed by evaluating Equation 14 , Equation 18 , and Equation 19 . This method will take computationally longer but it will calculate the collective motion of mobile species . It does however introduce two complicating terms, - 1 and c . - 1 will be addressed further but 39 determining the proper value of c , the concentration of is not clear. It is common to take the concentration to equal to the total number density of the mobile species, but this could be problematic if not all species are moving during the correlation time. It is possible that - 1 may also account for this concentration discrepancy , but more work needs to be performed to validate this assumption. The third method for calculating the conductivity would be to directly evaluate t he correlation with the charge of the species as shown in Table 4 . Now the concentration should be accounted for by direct computation of charge flux . However, this will be the most computationally expensive method , requiring sufficient simulation and correlation times. Ongoing work is evaluating the Q dependence of conductivity using this method in order to ensure convergent behavior is observed. At the completion o f this work, the best method for conductivity calculation is not clear, and will continue to be an area of interest for further research. 3 . 4 Dynamic s tructure f unctions : G(r,t), S (Q, ), & I (Q, t ) T he dynamics of a multi - component material with n different chemical species can be described by the total number density function, Equation 21 where r is an arbitrary position in the material, is the position of kth atom of species i at time t, and N i is the total number of atoms in species i. The correlation of the number density is the total van Hove correlation function 62 . ( Equation 22 ) 40 Equation 22 Equation 22 - part and - part Equation 23 The self - part calculates the probability of finding an atom after some time at a distance of . While the distinct part shows the probability of finding an atom when centered about a different species of atom. These functions show the dynamic relationship of lithium diffusion by combining the spatial relationship of pair - distribution function with th e motion of the mean squared displacement , allowing you to see a time evolution of the probability density . ii reports correlation and can be written as Equation 24 with g ii is N i /V . It is convenient to divide G ii by in order to obtain a normalization as is done in this work. The MSD is related to the von Hove correlation by , Equation 25 41 further, and can be approximated for oscillation - type dynamics by , Equation 26 where is the average displacement parameter . If the process is a continuous random walk then Equation 27 with D being the diffusivity of the atom. While the van Hove correlation function gives an intuitive perspective on how to imagine the probability of particle motion, real world experiments typically measure the features of the related by the four - dimensional Laplace transform of the G(r,t). Equation 28 function, tells us the intensity at which Q quasi - elastic neutron scattering (QENS) , on the BASIS instrument at Oak Ridge National Lab. In a QENS experiment the numbe r of scattered neutrons is measured as a function of the momentum energy transfer for collision events at specific scattering angle s . T he amount of energy transfer is in the range of one eV , resulting in features that manifest as a broadening of elastic peak . In our garnet system, incoming neutrons will scatter off mobile lithium coherently and incoherently relating to the self and distinct van Hove correlations . 119 42 Equation 29 the resulting spectrum will comprise of three major components. The elastic response comprising most of the scattered information, is modeled by a delta function where there is no change i n neutron energy after scattering off an element. A background function , which is applied to account for instrumental noise . Lastly, a Lorentzian peak , used to capture the small change s in energy transfer as a function Q and temperature . Taking the half wi dth half max ( H WHM (Q) ) of the Lorentzian, and plotting versus Q 2 , one can extract out the diffusivity based upon the shape of the plot, signifying different modes of diffusion. We looked at four different models as to describe the specific types of motion that could occur in the garnet sample. The most basic model would be that of continuous motion as 62 . At t he small Q limit where we are probing long time and distance scales, all models should reduce to the linear form of . If the motion of particles follows a jump diffusion model, as Q 0 applying a Fickian model is applicable when extrapolating fr There are three jump diffusion models investigated in this work. A Chaudy - Elliot model 121 (CE) , which assumes a uniform jump frequency and distance. A Hall - Ross model 122 (HR) , where the jump frequency is constant, and the jumps follow a gaussian distribution , designed to model diffusion within a confined volume. And a Singwi and Sjolander model 123 (SS), meant to capture the diffusion when a particle will alternate between vibrating about a 43 position, and directed motion. Equation 30 119 lays out the expressions used to model the three jump mechanisms. Equation 30 From these expressions it is possible to extract the residence time and is the mean jump length. We now have a method for directly observing the spatial jump length and residence time, features that independently are un obtainable through other methods. Each of the four models were fit to the QENS experimental data for L5LT at temperatures ranging from 50 K to 700 K. Figure 8 shows the fitting for L5LT to the four models Figure 8 : QENS HWHM fitting of L5LT at 600 K for possible diffusion mode ls. 44 at 600 K. Because of large coherent scattering at Q values larger than 1.2 Å - 1 our fitting range for termine the best model for capturing the jump mechanics between HR and SS. Clearly CE with its sinusoidal decay at large Q does not adequately apply, implying that the diffusive jumps are irregular in time and length. HR and SS models both predict a platea uing maximum at large Q, though our data does not allow us to easily see this feature. The major difference between HR and SS models relates to the definition of the average jump length distribution. HR models a gaussian jump distribution , in contrast SS assumes a gamma distribution and captur es more short hop ping events . A plot of vs Q 2 for the two models shows a more rapidly converging plateau for the HR model , while the SS model will gradually approach the maximum. Despite b oth models be ing reasonably fit to our data, the descriptive definition of the SS model more closely matches what we see in the garnet system and is the jump model used for determining the average jump distance and site residence time. A convenient means of relating G(r,t) and S(Q, function to make what is known as the intermediate scattering function I(Q,t). This partial transformation of G(r,t) by Equation 31 , Equation 31 t ranslates the real space component r to the Fourier space Q . Taking the Laplace transform of on the other hand changes the time variable from frequency to real time by Equation 32 . Equation 32 45 When modeling using molecular dynamics, it is possible to calculate t he species - specific incoherent intermediate scatteri ng function as , Equation 33 where is the number of atoms in species and is the position of nth atom at time t . The allowed values are Equation 34 where is the lattice parameter, the integers, and the unit Cartesian vectors. The neutron - weighted coherent ISF is than calculated according to Equation 35 where is the coherent scattering length (0 for shells) and is the total number of particles. The to tal neutron - weighted ISF is then the sum of coherent and incoherent contribution s, Equation 36 (arbitrarily weighted by the total incoherent scattering power), where is the number of species and the incoherent scattering cross section (0 for shells). The coherent and incoherent 46 scattering parameters of the elements based on the natural abundance of isotopes (except for Li where values from 7 Li isotope were used) are shown in Table 4 Table 5 : Coherent and incoherent scattering powers of elements involved in the present study. 7 Li La Ta Zr O (fm) - 2.22 8.24 6.91 7.16 5.803 (fm 2 ) 0.0078 0.0113 0.0001 0.0 00 2 0.000008 Because the signal from a QENS experiment takes the form of a Lorentzian peak for mobile species , taking the Laplace transform of the peak yields an exponentially decaying function in I(Q,t) . Fitting the exponential is equivalent to the Lorentzian fitting can be used as a means of calculating the self - diffusivity of the mobile species by evaluating the Q dependence of the decay . The relaxation of I(Q,t) for the garnet system did not appear to conform to strictly exponential decay and the stretched exponential model of Kohlrausch - Williams - Watts (KWW) 124 was applied . Equation 37 KWW is the relaxation f requency , and is the stretching parameter. Taking 1/ KWW will yield the relaxation time and the mean relaxation rate can be obtained from when G is the gamma function. and simulated I(Q,t) we re performed by from units eV to frequency and fitting th e HWHM from experiments to the calculated KWW . 47 3 . 5 Thermodynamic f actor O ne of the fundamental problems when running simulations is the relatively small number of part icles being simulated when compared to the number of species in an experiment. This difference of scale means it is likely that a property derived from the finite sampling using statistical mechanics may differ from one observed at the thermodynamic limit . 125 - 127 Such a property as discussed previously is the ionic conductivity . If a system is dilute and noninteracting, the inverse thermodynamic factor equals one, ( ) evaluating the conductivity of the system is equivalent to the basic Nernst - Einstein equation. As the particles in the system interact more and more, the valu e of will decrease proportionally to the degree of interaction. 128 When performing simulations on the garnet series , as done in this work, if you simply derive the diffusivity from the mean squared displacement and apply the Nernst - Einstein equation, the calculated conductivity will be overestimated by about an order of magnitude. is defined as the change in number density of particles with respect to a change in chemical potential. Numerically one can evaluate this expression as done in Equation 38 Equation 38 by selecting a fixed volume within a microcanonical (NVE) simulation and measuring the difference between the number of particles at some point in time with the average number of particles in the sub - volume. Effectively this sub - volume now models a grand c anonical system freely exchanging particles with those outside of the sub - volume. This results in energetically discrete configurations while the entire simulation remains at constant energy. Selection of an appropriate sub - volume was explored with the goa l of finding a sub - volume with enough species 48 to adequately provided minimal variations in the total number of particles, while still being small enough as to not envelope the entire simulation cell and provide no reservoir for particle exchange or exhibit self - reinforcing through the periodic boundary. To test the validity of an appropriate sub - volume, classical MD simulations were performed for 1 ns on L7LZ with varying simulation sizes ranging from a 1 x 1 x 1 unit cell, to a 4 x 4 x 4 supercell. A MATLAB script was generated in order to divide each simulation cell into smaller sub regions and track the particle count over th e last 500 ps of a simulation. As can be seen in Figure 9 the magnitude of - 1 is dependent on both the absolute size of the simulation as well as the size of the sub - volume within the simulation. As the size of the simulation increases it is possible to probe larger sub - volumes and thus provide a more appropriate representation of continuum behavior. In such an experiment, we hope to find converging behavior in order t o maximize the efficiency of the simulation and produce meaningful insights. Based upon other Figure 9 : Finite size effects of - 1 for L7LZ at 1400 K for 1 x 1 x 1 to 4 x 4 x 4 simulation cells. ( a ) Absolute magnitude of - 1 as a function of the inverse length of the sub - volume, solid lines denote guides to the eye. ( b ) - 1 as a function of the volume ratio of the sub - volume and the simulation cell with a third order polynomial fit. (a) (b ) 49 factors related to proper phase behavior, a 3 x 3 x 3 cell appears to have mostly identical properties to the 4 x 4 x 4 simulation cell and deemed sufficient for further evaluation . But as can be seen in Figure 9 there is still the question of what an appropriate sub - volume is even within a specific simulation size. For 1 x 1 x 1 the simulation cell is too small to reasonably provide statistical significance and this method should not be used. 2 x 2 x 2 simulation s allow for the inclusion of a 1 x 1 x 1 cell centered within the 2 x 2 x 2 total volume, and this appears to be a decent benchmark for such a sized simulation. In our first implementation of this technique, 2 x 2 x 2 simulations were sub divided into two sizes, one at 1/8 the total volume and 1/16 the total volume. These two values were averaged and applied to Equation 8 to find the final conductivity. A refinement o n this technique was attempted by increasing the number of sub - volumes within the simulation, until the entire simulation space was covered, averaging among the sub - volumes, and performing this for many possible sub - volume s izes . Figure 9 (b) , plots - 1 versus the ratio of the sub - volume to the total volume of the simulation cell. To account for the uncertainly in selecting a specific sub - volume, we attempted to fit - 1 with lines, third order polynomials, and fourth order polynomials in order to extrapolate - 1 at V sub /V cell = 0. T he re was also an attempt at linear fitting by select ing V sub /V cell values in the mainly linear region between 0.2 and 0.8, but ultimately f elt that the manual selection of an appropriate region was not rigorous enough to confidently predict conductivity values. For the third and fourth order polynomials, the entire V sub /V cell range was fit, however because of the rapid increase of - 1 at low V sub /V cell , it is likely that the extrapolated value of - 1 would under correct when applied to Equation 19 resulting in an over estimation of conductivity. 50 Kirkwood a nd Buff 129 first showed that thermodynamic properties such as partial molar volumes and activity coefficients can be solved for by integrating pair correlation functions in real space. Kr ger et al 130 expanded upon this by attempting to apply the theories of Kirkwood and Buff to liquid finite sys tems like those in a molecular dynamics simulation. The number density fluctuation can be obtained b y integrating g(r), the pair distribution function , with an analytic expression for the Dirac delta functions in three dimensions ( Equation 39 ). 130 - 131 Equation 39 In the case of the garnet system, with the assumption that the Li distribution is random , at large values of r we would expect g(r) to converge to one like a liquid solution. We attempted to extrapolate the value of - 1 as r by applying a window function and performing linear regression on the oscillating signal. As can be seen in Figure 10 (A) the g(r) for Li does not nicely converge to o ne within the range of simulation cells tested. It appears that the distribution of lithium has some form of order within the volume explored in this work . The r esulting plot shows nonuniform bumps in g(r) similar to that of a semi - crystalline material at large values of r. F ailing to converge neatly to one results in the decaying oscillatory signal obtained in Figure 10 (B) making the extrapolation at r difficult to determine. This technique while more rigorous than the fluctuation describe previously, does not appear applicable to the garnet as tested because the Li distribution is too dissim ilar to that of a liquid in which this approach was formulated. 51 Be ing limited to half the simulation size, makes the extrapolation as r goes to infinity quite difficult because there will always be a question as to whether making the simulation size any larger will result in more useful information at the expense of expedie ncy. One possible method of alleviating this problem is to obtain - 1 by taking the Fourier transform of g(r) or directly evaluating the static structure factor S(Q) at the limit as Q 0 ( Equation 40 ). Equation 40 This method allows for the use of the smallest Q vector possible for m a given simulation cell. Which in turn is longer than the half of a simulation cell we are limited to in the previous method and used to calculate - 1 for our core - only models. Figure 10 : Simulated Li Li PDF for L7LZ at 1400 K for various simulation cell sized. (a) (b ) 52 3 . 6 Excess e ntropy A s previously stated, it is thought that the transformation from tetragonal to cubic for L7LZ at high temperatures is driven by an increase in the entropy associated with going from an ordered to disordered lithium sublattice. Evidence to this point is inf erred from scattering experiments that show the lithium symmetry breaking about tetrahedral and octahedral sites but quantifying the change in entropy and determining whether it is a first or higher order phase transformation process is still debated. By a dopting the ideas of both Gibbs entropy and an atom being located in a specific region it is possible to compare the relative amount of disorder being introduced by the addition of heat or through substitution. 132 Information theory defines entropy as the measure of how much more information is required to define the exact state of a system 132 . For a statistical mechanical system this translates to determining the probability for a particle to occupy a microstate at that energy , the more particles and possible micro states , the more information required to define the system When performing a c onstant number, volume, and energy ( NVE ) MD simulation we are probing the configurational space at a spec ific energy , identifying possible configurational microstates. I f we track the position of each atom during the simulation, we can determine the probably that a particle will occupy some small volume called a voxel. Voxels are three dimensional pixel s and act as an analog to a microstate as only one atom can occupy a voxel at a time . I n this work each voxel is 0.1 Å in length , ensuring a constant vo xel volume at the expense of constant voxel number. 53 The probability of a voxel being occupied is calculated by summing the number of times a voxel is occupied divid ed by the total number of frames in the simulation. Plugging the for the configurational entropy. ( Equation 41 ) Equation 41 we then determine the entropy of a given garnet composit ion. Because each composition of garnet has different number of total lithium , absolute unit cell length, and voxel number making a direct comparison of entropy in each simulation cannot be performed. T o make a relative comparison possible we calculate am ount of normalized excess entropy generated by the addition of additional lithium. The normalized excess entropy is defined as Equation 42 where S ex is the normalized excess entropy, and S i is the ideal entropy. The ideal entropy is taken classically as the Boltzman entropy, whe re each voxel has equal probability of being occupied . A helpful analogy to explain how Shannon entropy works is the case of a fair and a biased coin. Plugging a fai r coin in to Equation 41 yields a result of one. While in the case of a coin with a 33% chance of being heads Equation 41 would be . In this system the amount of excess entropy would be - 0.18, meaning under ideal circumstances, a fair coin introduces more disorder than a biased one. In the limiting case of a double - sided coin S c = 0, and the state can be determined absolutely. We use the same type of evaluation in this work. The more ordered the system, the smaller the normalized excess entropy. As more lithium is added or when the temperatu re is increased, it is expected that the amount of entropy will increase . It is suspected that a 54 composition with a l arger excess entropy is correlated with a larger diffusivity and it may be possible to use excess entropy calculations in the future as a m eans of identifying highly conductive materials without having to perform the more computationally expensive correlation evaluations. We make a visualization of voxel probability during this process by generating a nuclear density map. The resulting file can be used to generate isosurface plots which represent the three - dimensional probability distribution function. These isosurface plots are like long exposure photographs allowing us to see the diffusion network for lithium by superimposing the positions occupied by the species throughout the simulation . By comparing the isosurface plots across temperature and composition for the garnet series we can visually see the changes in diffusivity and site occupancy allow ing for an intuitive look at the dynamic processes. 55 3.7 Statement on s tructure of this thesis W ith the background and methods established in Chapter 1 - 3 , Chapter 4 turns to the application of the acquired information to the lithium garnet series of ionic cond uctors. I begin by first looking into the phase transformation behavior of LLZ through DFT - MD simulations. It is critical that a model accurately captures the proper dynamic response upon heating in terms of lattice expansion coefficient and phase transfor mation behavior. As stated above the XC functional is a selectable parameter when running these types of simulations and choosing the right one for a system is typically a matter of trial and error. We took a look at 14 different XC functionals in this wor k, and while by no means encompass the total number of functionals in the literature, it made a representative sample for the most commonly used by other researchers. Because DFT is a higher - level simulation experiment, it is often used as the premier simu lation method for materials. Results from DFT are commonly used to form model parameters for lower level simulations like classical MD as done in this work. From this work we can gain a context as to how our lower level perform and evaluate sources of devi ation in our classical MD models. In Chapters 5 & 6 we evaluate the performance of the first classical MD model used in - dynamics involved in lithium motion. Chapter 5 compares and contrasts the behavior of the two end member compositions in Li 7 - x La 3 Zr 2 - x Ta x O 12 (x=0, 2) by first evaluating our ability to capture structural and diffusive behavior consistent with experimental results and then by evaluating the local environment of lithium, and how the differences in environment effect the macroscopic properties. T he compositions chosen represent two model materials within the garnet series, each containing only one type of element in the Ta/Zr site so combinatorial effects are negated. L5LT is the ideal cubic garnet, with lithium fully disordered about its sublatti ce, 56 contrasting with L7LZ with its ordered lithium sublattice and low temperature tetragonal phase. It is here we also introduce our first correlation function, the von Hove correlation function, laying the groundwork for the evaluation of diffusivity by Q ENS and the calculation of the intermediate scattering function. The properties of L5LT are expanded upon in Chapter 6 as we investigate the self - diffusivity through experimental QENS measurements and the simulated intermediate scattering functions. This novel work establishes QENS as a viable technique in the study of lithium - ionic conductors. The results of this study are meant to support the dynamic processes previously reported for the garnet series. Chapter 7 takes a look at a new simulation model, t - model was created to act as a more efficient means of calculating simulation trajectories by eliminating the calculations involving the shell in the previous work. This allows us to lengthen the timestep for each frame i n our simulation by no longer requiring us to capture the very fast electron polarization while simultaneously allowing us to use a more efficient simulation program utilizing a domain decomposition algorithm in way of a replicated data strategy. The valid ity of this new model is then tested on the model material L7LZ at various simulation sizes to confirm convergent property behavior. Lastly, Chapter 8 culminates our efforts in model selection and analysis techniques by looking at a range of possible lith ium concentrations within the garnet series. This work was performed in order to investigate the role of lithium concentration on the properties of the garnet in order to better explain the optimal composition for use in a solid - state battery. The previous core - shell model used in chapter 5 and 6 was only used to model two model garnet compositions. 57 Here the efficiency of the core - only model allows for seven different Li 7 - x La 3 Zr 2 - x Ta x O 12 (x=0, 0.5, 1, 1.25, 1.5, 1.75, & 2) compositions. 58 CHAPTER 4 : Effec ts of electron exchange - correlation functional s on the density functional theory simulation of phase transformation of fast - ion conductors: A case study in the Li garnet oxide Li 7 La 3 Zr 2 O 12 4.1 Abstract T he phase transformation of Li 7 La 3 Zr 2 O 12 (L7LZ) is investigated upon heating using first - principle molecular dynamics simulation by applying the local density approximation (LDA) and th irteen generalized gradient approximations (GGA) of the electron exchange and correlation (XC) energy functiona ls within the density - functional theory (DFT) framework. It was found that some functionals in the selected group failed to predict the phase transformation behavior while others predicted lattice volumes and lattice parameters larger or smaller than those experimentally. Of the fourteen, three functional types, PBEsol, SOGGA, and PBE2 , exhibited behaviors consistent with the tetragonal to cubic phase transformation upon heating and they were able to reproduce crystallite volumes with 1.5% of the experiment al values. The correlation of XC functional forms and their accuracy in predicting materials properties is discussed. 59 4 .2 . Results and discussion 4 . 2 .1 . Phase behavior at 1000 K: all 14 functionals T o first determine whether an XC functional merited further study, an initial structure of tetragonal L7LZ based on experimental neutron diffraction data 37 was heated to 1000 K to check if a tetragonal to cubic phase transformation could be observed, as the experimental transition temperature is aro und 900 K 97 - 98 . Figure 11 presents the results of these initial simulations by plotting lattice parameters (a, b, c) as a function of time. In general, lattice variation could be Figure 11 : Time dependent lattice parameter for the 14 tested functionals at 1000 K. 60 characterized as two categories. First, systems with functionals such as LDA, PBEsol - D3(BJ), VMT84, WC06, PBE - D3(BJ), AM05 and PBE - FE almost always stay as a tetragonal phase in the duration of simulation. These functionals are categorized as the tetragonal (T) group. The remaining functionals exhibit moment s of cubic phase before reverting back to some form of tetragonal structure. They are categorized as the T/C* group with varying degree of cubic/tetragonal volatility. Table 6 : Summary of simulation results using different XC functionals XC Functional Lattice Shape Lattice Volume Error (1000 K) 1000 K 900 K 700 K 300 K LDA 84 T - 2.86 - 3.53 - 3.51 - 3.57 PBEsol - D3 ( BJ ) T - 2.56 VMT84 104 T 0.40 WC06 105 T 0.46 PBE - D3 (BJ) T 0.77 AM05 75 T 1.46 1.50 1.30 0.84 PBE - FE 102 T 1.94 1.87 1.78 1.43 PBE3 100 T/C* - 1.02 - 1.25 - 1.00 - 0.01 SOGGA 101 T/C* - 0.36 - 0.43 - 0.45 0.72 PBEsol 105 T/C* 0.31 0.35 0.23 - 0.01 RGE2 116 T/C* 2.70 PW91 113 T/C* 3.87 3.81 3.49 2.97 PBE 86 T/C* 4.30 4.12 3.87 3.35 RPBE 115 T/C* 8.17 61 In Table 6 , we present the quantitative comparison of lattice volume at 1000 K obtained from simulation to that of extrapolated X - ray diffraction data from Matsuda 98 . We can see that lattice volumes are delimited by LDA (smallest) and RPBE (largest), which can be correlated to h ow the exchange energy depends on the reduced gradient, to be shown in Figure 16 . LDA and PBEsol - D3(BJ) functional underestimate the volume by more than 2 %, while PBE, PW91 and RPBE functionals overestimate the volume by more than 3%. Some of the newer functionals specifically designed for solids, e.g. AM05, WC06, SOGGA, PBEsol are able to predict volume values within 1.5%. These two observations were also made in several other studies of effect of XC functionals 89 - 96 . Other more recent functionals such as VMT84, PBE - FE, PBE2, and RGE2 have scattered success in obtaining values consistent with experiments. While VMT84 predicts a very accurate lattice volume, it does not predict the right phase behavior. Both PBE - FE and R GE2 slightly overestimate while PBE2 slightly underestimates the volume. When we compare dispersion corrected and non - corrected functionals, it seems that addition of dispersion energy reduces the lattice volume, which brings lattice volume from PBE - D3(BJ) close to experimental values and leads to underestimation from PBEsol - D3(BJ). This observation is similar to those reported in the literature on strongly and weakly bound materials 96 . 4 .2 .2 . Lattice parameters at a ll temperatures: 8 functionals B ased on results in Figure 11 / Table 6 and popularity of functionals in the field, we further selected 8 functionals, i.e. AM05, LDA, PW91, PBE - FE, SOGGA, PBEsol, PBE, and PBE2, to perform simulations at lower temperatures (300, 700, and 900 K) to investigate the phase transformation and latti ce volumes. Plots of the time dependent lattice parameters at 300 K, 700 K, and 900 K, can be found in Figure 13 , Figure 14 , and Figure 12 respectively. 62 The 900 K trials are of particular interest because the response of the cell at or near the transition temperature can very widely throughout the time of the simulation. It appears the eight functionals that we decided to evaluate further remain consistent with their 1000 K analogs with a major exception being the response of LDA. In the 1000 K simulation LDA consistently maintains a tetragonal structure over the length of our averaging time, however it is clear at 900 K LDA is capable of capturing the phase transformation, its susceptibility to over - b i nding st ill makes it a poor candidate as an XC functional to use in the garnets. The 300 and 700 K simulations well below the transition temperature for this material show consistent tetragonal phases in agreement with our expectations. It is evident from comparing the lattice parameters at these temperatures which models will over - bind (LDA), and which will under - bind ( AM05, PBE). The other five functionals evaluated overall have quite good agreement with experimental measurements. Figure 12 : Time dependent lattice parameters for 8 XC functionals at 900 K 63 Figure 13 : Time dependen t lattice parameter for 8 XC functionals at 300 K Figure 14 : Time dependent lattice parameters for 8 XC functionals at 700 K 64 Figure 15 consists of the time averaged lattice parameters with standard deviation for L7LZ at all four simulated temperatures. Each simulation is then compared to the X - ray diffraction structure reported by Matsuda et al. 98 . Their study shows a transformation temperature of about 927 K, consistent with our experimental investigation of impurity free L7LZ 97 . While all structures in Figure 15 maintain consistent tetragonal character below 900 K, we consider that four functionals, i.e. SOGGA, PBEsol, PBE, and PBE2 functionals are able to predict a cubic - like average structu re at 1000 K due to the overlapping standard deviations. Out of these four, SOGGA, PBEsol, and PBE2 are able to predict lattice volume within 1.5% of experimental values at all temperatures. 65 Figure 15 : Average lattice parameter s for 8 XC functionals at 300 K, 700 K, 900K, and 1000 K (red). Experimental lattice parameter by X - ray diffraction by Matsuda. 66 4 . 2. 3 Correlation between XC functionals and phase behavior/volume I t is common in the literature to plot the density gradient dependence of the exchange enhancement factor to examine the difference of XC functionals, since there is generally less variation in the correlation energy. The plot of different GGA functionals employed in this study is shown in Figure 16 . The second - order gradient expansion (GE) of electron gas with slowly - varying density is shown as the dotted line. In general, it is proposed that a functional will overestimate the volume if its enhancement factor is larger than the gradient expansion at small s. This is the case for RPBE, PBE, PW91, and PBEFE. On the other hand, functionals with similar values to the dotted line at small s ge nerally can predict volumes close to experiments, e.g. PBEsol, WC06, SOGGA, PBE2, with the exception of RGE2 and VMT84, both of which have higher values than the other four at large s. This suggest that functional forms at both small and large s range affe ct the predicted lattice volumes. This is also illustrated with the AM05 functional whose deviation from the GE line at small s is likely to be compensated by its large value at large s. Out of the four functionals with good accuracy in lattice volume, PBE sol, Figure 16 : Plot of exchange enhancements factors of GGA functionals employed in the study. The second order gradient expansion as is shown in the dotted line. 67 SOGGA, PBE2, and WC06, the first three are also able to predict the right phase behavior, possibly due to their lower values at large s compared with WC06. We think part of the difficulty of these GGA functionals to predict correctly both the phase tr ansformation behavior and lattice volume is the complex bonding environment (a mix of ionic and covalent bonding with different elements) in Li 7 La 3 Zr 2 O 12 , along with the intrinsic limitation of GGA approximations as they only depend on the gradient (1st or der) of density. It was proposed that approximations involving higher - order terms (Laplacian of density), i.e. meta - GGA, could lead to more accurate results within the DFT framework, as demonstrated with the Minnesota meta - GGA MN15 - L 133 and Strongl y Constrained and Appropriately Normed (SCAN) functionals 134 . We plan to investigate the effect of these meta - GGA functionals in the future. 4. 3 Conclusion T he present study shows that molecular dynamics modeling of phase transformation of fast - ion conductors within the DFT frame - work can be a challenging task as it strongly depends on the employed XC functional. For the PBE - like GGA functional , the functional forms at both small and large density gradient ranges affects whether the correct lattice volume or phase behavior can be obtained. Of the fourteen functionals studied, SOGGA, PBEsol, and PBE2 show strong agreement with experimental crysta l structures of L7LZ and appear to be candidates for further study into the dynamics of lithium diffusion in this material family. 68 C HAPTER 5 : Core - s hell M odeling of L5LT and L7LZ 5.1 Abstract T o better understand the processes that determine the ionic conductivity and phase transformation behavior of the lithium garnet series we employ molecular dynamic simulations on two model materials within the series L5LT and L7LZT. We first compare the la ttice parameter, and phase behavior using NPT simulations, and the conductive/diffusive properties using NVE type simulations to experimentally measured values . We find good agreement for structural and dynamic properties of both compositions when a thermo dynamically corrected Nernst - Einstein equation is applied . We show that a general conduction pathway is difficult to determine because trajectory of lithium through the shared triangular bottleneck conjoining the neighboring polyhedral is observed between the two compositions and changes with temperature. We also show that at higher temperatures and concentrations of lithium, while at low temperatures most lithium atom s perform an oscillatory jump diffusion mechanism. The differences in diffusion mechanism appear to be related to the local arrangement of lithium with its nearest neighbors, nter - type trajectories and the asymmetrical clusters of L5LT exhibit - trajectories through the shared bottleneck. 5.2 Results and discussion 5.2.1 Validation of simulation process T o demonstrate the validity of our - force - field model and our simulation approaches, we compared lattice parameters, neutron scattering PDF, and ionic conductivities calculated from MD simulation with those from literatures for both L5LT ( Figure 17 (a)) and 69 L7LZ ( Figure 17 (b)) . The L5LT model presents lattice parameters that match well with X - ray and neutron di raction data independently collected from Wang et al. 37 , Cussen 3 , Thangadurai et al. 25 The total scattering PDF ( Figure 17 (a)) from neutron diffraction collected by Wang et al. 37 also matches well with simulation results suggesting our simulations provide a representative structure of cubic garnet. The L7LZ core - shell model predicts lattice parameter and thermo - expansion coefficients smaller than those measured experimentally . 37, 39, 43 The thermal expansion along the a/b directions in particular is smaller than expected in the tetragonal phase, and after transformation to cubic, still grows at a smaller rate than measured by Matsuda. 98 The tetragonal phase appears to constrict above 700 K in our simulations contrasting the experimental measurements which show a consistent growth until the transformation temperature is achieved . The rate of ther mal expansion along the c direction increases in a nonlinear fashion until the temperature reaches 900 K, above which L7LZ is in the cubic phase. When the tetragonal lattice parameters were converted to pseudo - cubic ones , by taking the cube root of the vo lume , a linear expansion was observed for the whole temperature range. This transition temperature of 900 K agrees with the reported transition temperature of 918 K by Larraz et al. 43 , 9 13 K by Matsuda et al. 98 and 923 K by Matsui et al. 4 using high - temperature XRD and thermal analysis. First principles DFT work by Bernstein et al. 44 predicted a transition temperature between 800 1000 K but did not investigate the transition temperature to a higher resolution. This temperature of 900 K is much higher than those reported elsewhere, e.g. 373 423 K from Geiger et al., 30 450 K from Adams et al., 50 623 K from Kuhn et al. 31 I t is plausible 70 that these samples were contaminated with Al, H 2 O, or CO 2 , which lowers the temperature of phase transition as shown by Wang . 97 The comparison of PDF at 300 K, in the form of G PDF (r), from MD simulation and neutron scattering experiments is shown in Figure 17 ( d) . Since we do not have neutron total - scattering data for L7L Z , we calculated G PDF (r) based on the Rietveld refinement results of conventional neutron di raction data from our previous work. 37 Specifically, we utilized the Figure 17 : (a) Lattice parameter of L5LT compared to neutron and x - ray diffraction studies. (b) Atomic pair distribution f unction from total scattering experiments compared to MD simulations (a) (b ) ( c ) (d ) 71 software package PDFGUI 126 to calculate G PDF (r) from the structure file of L7LZ at 300 K, supplied as supplementary data in the literature. 37 approach. Because Rietveld refinement assumes independent atomic models, it is expected that this approach will generate PDF patterns similar to those from total - scattering experiments but with variance in peak inte nsities and width. Such variance is due to the di erence in the average structure (from the Rietveld refinement) and the local structure, as shown in the study of Bi 2 O 3 . As shown in Figure 17 (d) , the PDF plot from MD simulation agrees well with that from the suggesting atomic correlation. Finally, we looked at the ionic conductivit y values for L5LT and L7LZ as a function of temperature to ensure not only accurate structural properties but also dynamic ones . Figure 18 (a) & (b) shows the MD simulation results for L5LT and L7LZ respectively A pplying the Nernst - Einstein (N - E) relation directly (blue lines) predicts ionic conductivity values an order of magnitu de larger than those observed experimentally 4, 25, 35, 97 but will within the reported range of other simulations that did not perform a thermodynamic correction. 5, 47 The inverse thermodynamic factor - 1 for these compositions is shown in Figure 18 (c) as calculated from the fluctuation method ( Equation 36 ) . The addition of - 1 bring s the conductivity val ues to well within experimental reports particularly with respect to the high temperature measurements though the low temperature conductivities for both models are divergent . Matsui et al. 4 and Wang et al. 97 both report a superionic transition 135 for L7LZ , where a poorly conductive phase turns into a more conducting phase upon heating. The superionic transition temperature coincides neatly with the tetragonal to cubic phase transition in Figure 17 (b) at 900 K . 72 It is possible that experimental samples contain more disorder than our simulation cells and this contributes to the higher conductivities. Another classical MD simulation, labeled as but predicted the transition temperature to be 400 K (not shown). This simulation also yielded Li Li pair distribution function inconsistent with experiments, which will be discussed later. Figure 18 : (a) conductivity for L5LT and (b) L7LZ as a function of temperature. Blue lines represent the conversion of MSD to D*, and then applying the Nernst - Einstein (N - E) equation ( Equation 8 ). Red lines are the thermodynamically corrected ( - 1*N - E) Nernst - Einstein equation. (c) Comparison between the inverse thermodynamic factors ( - 1) for L5LT and L7LZ using the fluctuation method ( Equation 38 ). L5LT L7LZ (a) (b ) ( c ) 73 The failure to obtain accurate conductivity values at low temperat ure is a major area of exploration in ongoing research. The decrease in conductivity for L5LT reported by Thangadurai 26 , is mostly related to the properties of the total conductivity , incorporating grain boundary resistance, and the possible addition of impuritie s. However, Wang in his bulk conductivity measurement actually reports a slightly higher conductivity than that obtained from simulation. For L7LZ it is much of the same story, the activation energy extracted for m the simulations appear much more in line w ith the first few data points, but there is clearly some form of nonlinear change in ac tivation energy as the temperature continues to decrease. At present we do not have a means of rectifying the discrepancies between the simulations and experimental values, though investigations leading us to peruse more sophisticated experimental and analysis techniques are expanded upon later in this thesis. The low values of - 1 observed in our simulations indicate each Li hopping e vent is not independent and is a ected by how the other Li in the system are moving. The inverse thermodynamic factor for L7LZ , reported in Figure 18 is consistently low er than that calculated for L5LT . 51 This is a reasonable result considering that increasing the Li concentration limits the number of ways in which Li will locally distribute themselves. This point will be discussed further in the later section regarding lithium clusters and dynamics. It ca n be seen from Figure 17 and Figure 18 that our MD simulations have predicted physic al parameters that agree reasonably with those from experiments. With this conformation , i n the subsequent sections we will take a closer look at the atomic - scale details of local structure and dynamics for L5LT and the phase transition of L7LZ . 74 5.2.2 L7LZ van Hove c orrelation f unctions (VHCF) Figure 19 shows (a) the self - part of VHCF as 4pr 2 G s (r,t) of lithium atoms and (b) its projection along the space and time axis at 300 K. There is a high probability of finding lithium atoms with displacement of 0.2 Å and such probability shows weak time dependence for di erent displacement values (0.11, 0.23, 0.58 Å). This type of correlati on is typical of atomic oscillation around equilibrium positions, which is expected for lithium atoms in the ordered tetragonal phase. Figure 19 also shows the (c) distin ct - part of the VHCF as G d (r,t) of lithium atoms and (d) its projection along the space and time axis at Figure 19 : (a) Self - part of the van hove correlation function and (b) its spatial and temporal projection at 300 K. (c) Distinct - part of the van Hove c orrelation function and (d) its spatial and temporal projection at 300 K. 75 300 K. The dependence on the distance is characteristic of pairs showing short - range order (peaks at small distances) but long - range disorder (oscillati ng around 1 at large distances). Such structure features are generally found in liquids or glasses. A small time - dependence at short time, e.g. 0.1 ps, was observed. The high - level presentation of VHCF is shown in Figure 20 as MSD and PDF plots for both 300 K and 1100 K. Again, a flat MSD plot at 300 K in Figure 20 (a) is a signat ure of oscillation. The linear MSD plot at 1100 K is a signature of di usion that we will discuss more in Figure 21 . In Figure 20 (b) , we also present the Li literature. Using the same approach discussed previously for the neutron scattering PDF, we Li PDF and we found it is similar to that from our MD simulation. Again, the slight difference is due to the difference between the average and local structure, as discussed previously. However, MD simulation by Adams et al. generates two peaks (2.35 and 2.8 Å), sug gesting that their Morse - type force - fields may not be able to capture the real (a) (b ) Figure 20 : (a) Mean squared displacement (MSD) for LLZ at 300 K and 1100 K. (b) Partial pair distribution function of Li Li pairs at different temperatures compared to that of Adams et al. 76 dynamics of lithium motion. With the increase of temperature from 300 to 1100 K, the first peak in the Li Li pair becomes wider and moves to a smaller distance, indicating weake r Li Li repulsion at higher temperatures. Figure 21 shows (a) the self - part of VHCF as log 10 4pr 2 G s (r,t) of lithium atoms and (b) projectio n of 4pr 2 G s (r,t) along the space and time axis at 1100 K. It is apparent that lithium atoms at this temperature can move di erent distances for di erent time intervals, as expected for a di usion mechanism. To understand such complex time and space depende nce and di usion details, we can look at a simple model system. A linear MSD plot at 1100 K as seen in Figure 20 (a) is consistent with a continuous random - walk or structureless di usion model. Based on the di usivity extracted from Figure 20 (a) , we then plot the self - part of VHCF of a continuous Figure 21 : (a) Self - part of the van Hove correlation function and (b) its spatial and temporal projection at 1100 K. (c) Self - part of the van Hove correlation function and (d) its temporal projection of a continuous random - walk model at 1100 K. 77 random - walk model in Figure 21 c and d. As expected, there is one peak along the distance ax is in Figure 21 (d) due to the competing quadratic and exponential term in Equation 26 . Similarly, there is one peak along the time axis due to the competing power and exponential term in Equation 27 ). In Figure 21 (c) compare Figure 21 (b) and ( d ) , it is clear that selected projection plots at di erent spatial and temporal values in Figure 21 (b) deviate from this continuous random - walk model. For example, - in the overall VHCF plot in Figure 21 . When we compare Figure 19 (a) and Figure 21 ( a ) and ( c ) , we can propose that di usion behaviors at 1100 K, i.e. Figure 21 (a) , are intermediate between an oscillation type and continuous random - di ielastic neutron scattering (QENS) experiments were used to study atomic di usion in materials. Several models, e.g. Chudley Elliott, 27 Hall Ross, 28 Singwi Sj ö lander, 29 has been proposed to quantify the QENS results. Parameters in these models are related to VHCF parameters by taking the Fourier transform as examined in Chapter 3 . A Look into these jump models is expanded upon in the Chapter 6 of this the s is. 78 5.2.3 Lithium distribution and its relation to phase transition Probability density functions (p.d.f.) or density maps o e r an intuitive way to visualize lithium distribution in the tetragonal and cubic phase. Figure 22 (a) shows the isosurface plot (isosurface level of 0.3 Å 3 ) for a slice a long the (111) lattice plane that constitute the looping structure of Li sites in the Li di usion pathway, where the two tetrahedral sites (Td - 8a and 16e) and two octahedral sites (Oh - 32g and 16f) are all visible. Figure 22 ( a ) is based on Rietveld refinement results so the density maps are local and harmonic (spheres or ellipsoids). Our 300 K MD simulation ( Figure 22 (b) ) predicts lithium densi ty to be primarily located at the Td - 8a, Oh - 32g, and Oh - 16f sites in agreement with Figure 22 (a) . However, we also observed a degree of anharmonicity arising from an elon gation of the Oh - 16f and Oh - 32g sites toward the Td - 16e sites. This elongation of the Oh - 32g sites that connect Td - 8a to Td - 16e continues to grow as the temperature is increased, leading to a small occupancy of the Td - 16e site at 500 K ( Figure 22 (c) ). At 800 K ( Figure 22 (d) ) a near complete loop structure is observed for the tetragonal Figure 22 : Lithium density maps on (111) for LLZ. (a) 300 K by Revitfield refinement, (b) 300 K, (c) 500 K, (d) 800 K, and (e) 1100 K by MD simulation. Isosurface level is 0.3 Å - 3 79 crystal. After the phase transition to cubic has occurred, the probability of occupying either of the now degenerate tetrahedral sites is equal, and a complete ring structure is observed as shown in the density map at 1100 K ( Figure 22 (e) ). Additional insight into the phase transition can be obtained by examining the cage occupancy and MSD of lithium at di erent sites shown in Figure 23 . The change in occupancies of tetrahedral sites with increasing tempera tures ( Figure 23 (a) ) clearly indicates the depletion of Td - 8a site and filling of Td - 16e, especially around the transition temperature of 900 K. On the other hand, the occupancies of the two octahedral sites remain close to one at all temperatures but do experience a decreas e in occupancy as the temperature is increased. Equivalently, there is a small increase in the overall occupancy of the tetrahedral sites. The shift of lithium atoms from Td - 8a to Td - 16e site was also proposed by Bernstein et al. 44 but they stated there was n o exchange of lithium atoms between tetrahedral and octahedral sites. In the meantime, slightly lower Td occupancy was observed in their work, as shown in Figure 23 (a) . Error bars in this Figure 23 : (a) Tetrahedral and octahedral cage occupancy for L7LZ as a function of temperature. Error bars represent the degree of dynamic fluctuation in the occupancy observed during MD simulation. Values from Bernstein et al are shown for comparison. (b) MSD for select ed atoms occupying three different cages type at 300 and 500 K. (a) (b ) 80 figure represent the degree of dynamic fluctuation in occupancy observed during MD simulation. Furthermore, the MSDs of all three sites at 300 K in Figure 23 (b) indicate lithium atoms oscillating as shown in the VHCF plot of Figure 21 . At 500 K, MSDs of Td - 8a and Oh - 32g lithium atoms show a simultaneous increase at short time suggesting that their motion is correlated. From the above observation and discussion of both Figure 22 and 23 , we believe that the tetragonal to cubic phase transition is an entropy - driven one that involves redistribution of lithium atoms among all tetrahedral sites. The ordered tetragonal structure is energetically favored as it has the smallest number of nearest - neighbor Li Li repulsion pairs. 37 Redistribution of Li among all tetrahedral sites will increase the entropy but is likely to raise the internal energy, since an occupied Oh site (3 2g or 16f) has a maximum of only two Li Li pairs while an occupied Td site (8a or 16e) has a maximum of four Li Li pairs. Above the transition temperature, the entropy contribution overshadows the energy contribution. The transition is likely to initiate o n the Td - 8a site but needs the cooperation of neighboring Oh - 32g sites. Lithium atoms at Oh - 32g sites mainly act as relay atoms, instead of contributing directly to the filling of Td - 16e site. Again, the reason that Oh sites prefer to be consistently occup ied is to minimize Li Li repulsion. However, there is indeed a slight shift of lithium atoms from octahedral to tetrahedral sites with increasing temperatures, since the higher temperature is more accommodating to increased Li Li repulsion. 81 L5LT with is lower lithium concentration and constant cubic phase does not exhibit the lithium sublattice ordering as descr ibed for L7LZ . Instead its lithium distributions are random and the density maps of L5LT show much more interconnectedness between the Td and Oh sites. Figure 24 shows the Li density at various temperatures for L5LT along the [100] direction as well as for slices along the (111) plane. Importantly from these projections we do not see any lithium density directly linking between Oh - Oh cages, instead lithium must pass through the Td site in order contribute to the conduction. This is consistent for both L5LT and L7LZ showing that there is no change in jumping mechanism between the two compositions . A ny difference in the magnitude of the conductivity must result from the degree of interconnected ness between lithium and their nearest occupied or unoccupied site. Figure 24 : Lithium density maps for L5LT derived from MD simulation at (a) 700 K, (b) 475 K, and (c) 300 K along the [100] direction. Isosurface level of 0.1 Å - 3.Two - dimensional Li density maps for the (111) plane, with a distance of 23 Å to the origin, cutting through Td - Oh - Td cages at (d) 700 K, (e) 475 K, and (f) 300 K with isosurface levels from 0 to 1 Å - 3. 82 5.2.4 Lithium clusters We will now take a closer look at the structure of nearest - neighbor Li Li pairs and their effects on lithium dynamics. There are theoretically 16 different ways in which lithium can distribute itself about neighboring Td and Oh sites, des ignated as clusters. A simplified notation has been previously introduced to easily differentiate the types of lithium clusters. 51 First, the is then used to identify a site as eith notation refers to the number of nearest occupied lithium sites. For Txx this ranges from 0 4 because there are four Oh sites surrounding a Td site, and for Oxx this ranges from 0 2 with two neighbo ring Td sites. Figure 25 : Possible local arrangements of Li with respect to their nearest neighbors 83 Figure 26 shows the evolution in the number of tetrahedrally and octahedrally centered clusters with temperature. It is expected that e ach cluster is associated with a unique potential energy and that the exact conduction path of a hopping event is decided by the ty pe of clusters that the hopping Li belongs to at the initial and final states. The number of clusters was evaluated Figure 26 : Average number of each type of possible lithium nearest neighbor arrangement as a function of temperature for L5LT centered upon Td (a) and Oh (b) sites, LLZ for Td (c) and Oh (d). ` L5LT L7LZ (a) (b ) ( c ) (d ) 84 from the last 500 ps of simulation, with t he error bars in Figure 26 represent ing the standard deviation of each cluster type at a particular temperature. Looking at L5LT first, at 300 K the most common cluster s for tetrahedral sites are the T11, T12 , T04 , T03, and T10 types . These particular arrangements suggest that the system minimizes the total Li Li interaction by maximizing the number of octahedral Li around empty an empty Td and minimizing the number of Li around an occupied Td cage. It is for this reason, the T14, an d T 00 have the lowest probability of occurrence. The two most common occupied Td centered sites, the T11 and T12 , can be seen as asymmetrical clusters that lead to uneven Li Li interactions that will offset both the Td and Oh lithium from their crystallog raphic centered position. This in turn will drive the dynamics of lithium diffusion to take what is called and - y where the lithium diffuses close to the bottle neck edge joining neighboring T d and Oh polyhedral. As the temperatur e is increased, we see the number of T11 and T04 configurations dramatically drop, while T03 and T0 2 distributions see a proportional increase. T he overall distribution of configurations becomes more uniform at higher temperatures, where the increase in kinetic energy is sufficient to overcome unfavorable interactions. Coupled with the distribution is the increase in standard deviation with each co nfiguration. The fluctuations of clusters is related to the number of hopping events, so the degree of fluctuation indirectly indicate the stability of each cluster type at a given temperature. At high temperatures the lithium is equally mobile across all observed cluster types as expressed from the large error bars. The analysis of cluster distribution has been focused on the T x x type only, though similar reasoning can easily be applied to the trends observed in the Oxx types. 85 Figure 27 ( a ) shows a sample of the different dynamic processes observed in our simulatio ns for L5LT . In the fist block a O11 cluster is shown with lithium rattling about their crystallographic sites. There is a very large amount of displacement observed for each site with the tetrahedral moving about almost the entire volume of the cage. The electrostatic repulsion between the neighboring lithium causes the Oh Li to be heavily displaced from its ideal site, only moving to the middle of the cage after the tetrahedral has displaced itself to the far side of its cage. There are also two types of single particle diffusion events observed in Figure 27 ( b ) and ( c ). The first type of event shown in Figure 27 ( b ) is of a Td to Oh hopping even. Here a lithium moves from a T12 or T13 type cluster into a neighboring unoccupied Oh site. This process is Figure 27 : Observed dynamical events at 500 K for L5LT. Yellow circles represent Li ions. Pink squares are tetrahedral sites and blue rectangles are octahedra lcages. In each example, the local environment, Li trajectories, and geometries of bottlenecks as Li goes through the faces are illustrated. (a) (b ) ( c ) (d ) 86 carried out by lithium diffusing through the bottle neck as shown in the accoupling diagram. This - - pass diffusion is shown in Figure 27 ( c ), when a single T03 type lithium diffuses by a Oh - Td - Oh diffusion pathway . It is clear in this third example how the clustering of lithium impacts the overall trajectory followed, the non - diffusing T03 lithium act as a repulsive force on the diffusing lithium forcing it to pass through the edge of bottleneck and immediately into a neighboring unoccupied This minimizes the Li Li interactions by avoiding the formation of an unfavorable T13 type cluster. The final type of diffusion event presented in Figure 27 d ) is an example of collective diffusion , were two neighboring lithium diffuse simultaneously through the center of their respective 136 AIMD calculations on the L7LZ composition and expected to be more common as the concentr ation of Li increases due to the interconnectedness of the sublattice and increased symmetry of the lithium clusters . The dynamics of clusters for L7LZ tell a somewhat different story. For the low - temperature tetragonal phase at 300 K, there are only two t ypes of observed clusters, T14 and T04 resulting from ordering about lithium sites, to minimize Li Li repulsion. As we increase the temperature, the number of T14 and T04 clusters decreases, first leading to an increase in the originally unobserved T13 and eventually the creation of T03 clusters. A similar trend can be identified for the Oh clusters. Initially only O11 and O10 clusters are observed. As the temperature increases the amount of O10 remains relatively constant, but the concentration of O11 decr eases significantly while we see the generation of the previously unobserved O12 and O02 clusters. It is to be noted that it was reported by Xie et al. 137 and Jalem et al. 47 that it was not possible to have O 12 clusters in L7LZ , based on the exclusion principle. Error bars in Figure 87 26 show the degree of fluctuation over the course of a simulation, with larger error bars signi fying larger fluctuations from more mobile ions. The dominance of clusters with high coordination numbers, e.g. Tx3, Tx4, implies that highly correlated motion is expected, consistent with the low inverse thermodynamic factors throughout the temperature ra nge, and with lower inverse thermodynamic factors in the tetragonal phase as seen in Figure 18 (c) . T o visualize the oscillatory and di usive motion discussed in Figure 19 and Figure 21 . Figure 28 ( a ) and ( b ) show the schematics and trajectories (connected small beads) of vibratory motions exhibited by lithium atoms in T14 and T04 clusters at 300 K. Within a T14 cluster ( Figure 28 (a) ), the center is Td - 8a site, and four neighbors are Oh - 32g sites. Lithium atoms at Figure 28 : Examples of lithium dynamics at 300 K for (a) T14 cluster, (b) a T04 cluster, and at 1100 K for a (c) T14 cluster, and (d) a T04 cluster. Yellow dots represent Li atoms. Pink squares and blue rectangles schematically represent Td and Oh cages 88 Oh - 32g site s are rattling around positions that are far away from the occupied Td - 8a. Wit hin a T04 cluster ( Figure 28 (b) ), the center is Td - 16e site, and two neighbors are Oh - 32g and the other two are Oh - 16f sites. Lithium atoms at Oh - 32g and Oh - 16f sites ar e rattling around but staying close to the empty 16e site. Both types of oscillation can be considered as ways to minimize Li Li repulsion. Figure 28 (c) and ( d ) show the schematics and bottlenecks of diffusion in L7LZ at 1100 K. In Figure 28 (c) , two lithium atoms (Li1 and Li2) will move in concert to change a T14 cluster to a T0 4 cluster. During the process, two triangular bottlenecks that connect the neighboring tetrahedral and octahedral cages were crossed. Locations of lithium atoms on these bottlenecks are also shown. The observation that lithium atoms prefer to stay close to one of the - suggested by the NEB study by Xu et al., 36 FPMD simulation by Jalem et al. 47 and Miara et al. 5 , and our classical MD simulation for Li5LT at low temperatures. 51 Figure 28 (d) presents the correlated motion of three lithium atoms crossing four bottlenecks. First, both Li1 and Li2 move synchronously with Li1 expelling Li2 into the empty Td cage. Effectively this changes a T04 cluster to a T14 cluster. Then Li2 and Li3 move synchronously in that Li3 migrates into the empty Td cage on the left while Li2 hops into the Oh cage that hosted Li3 previously, - - To obtain statistical detail s on how lithium atoms move across the bottlenecks, we recorded the Li locations on the bottleneck for a simulation time between 900 ps to 1 ns in the MD simulation. More than 1000 lithium atoms were collected we used both the distance to three corners, de noted as on - face - Li - to - oxygen distances, and distance to three edges, denoted as on - 89 face - Li - to - edge distances, to project where lithium is relativity located as it di uses through the bottleneck. Figure 29 compares the distribution of these distances between L5LT and L7LZ at various temperatures with a comparison to reverse Monte Carlo simulations performed by Wang 37 . For both compositions a Gaussian distribution was observed for t h e lithium to oxygen distance , however the lithium to edge distance esp ecially at low temperature s presents a bimodal distribution. As the temperature is increased, the bimodal distribution becomes Gaussian in both cases as the lithium clusters redistribute to obtain a more uniform distribution. L5LT interestingly exhibits more edge - pa ss type events signified by the increases fraction of events at low Li to edge distances. - - understood from the consideration of global and local symmetry. It is expected that t he centers of a tetrahedron and of a bottleneck would be favored for Li to pass from the perspective of global symmetry, while local symmetry might distort this picture, especially at lower temperatures. In Figure 26 (b) we saw that T04 dominates the T0x clusters in L7LZ . It can be reasoned that the local symmetry of T04 clusters causes a more uniform repulsive force on the Li, making the center of the bottl eneck the easiest for the Li to approach. In Li5LT there was a greater variety of clusters observed especially at low temperature. These clusters have a greater local asymmetry at most tetrahedrons due to Li + not fully occupying the surrounding octahedral sites, causing an asymmetrical force to be applied to Li and push it to the edge of the bottleneck. At higher temperatures, the influence of local Li Li repulsion becomes less important and Li will proceed - consistent with the global symmetry. 90 L7LZ Figure 29 : Histograms of lithium position on the face of the bottle neck with respect to distance from the nearest oxygen or shared edge of oxygen polyhedral for LLT and LLZ. (a) dia gram of bottle neck with solid lines representing distance to nearest oxygen, dashed lines representing distance to nearest face. Shown are events for LLT at (b) 550 K, (d) 900 K, and (f) 1100 K. Shown for LLZ at (c) 600 K, (e) 900 K, and (g) 1100 K. (a ) (b ) (c ) (d ) (e ) (f ) (g ) L5LT 91 5.3 Conclusions U sing molecular dynamics, we investigated local structure and dynamics in the solid electrolyte Li 7 La 3 Zr 2 O 12 and origin of tetragonal to cubic phase transition. First, validation of our force - field parameters was accomplished by comparing our simulation output to lattice parameter, neutron scattering, and conductivity measurements obtained from experimental tech niques. Thermal expansion, phase transition, and neutron scattering PDF, and conductivity results at various temperatures concur with experimental values. Second, we show that lithium dynamics predominantly consisted of oscillations at low temperatures in the tetragonal phase and combining self and distinct part of van Hove correlation function and MSD data to make a more complete analysis of structure and dynamics of Li. Third, we show that the tetragonal to cubic phase transition coincides with the redistribution of Td - 8a lithium to Td - 16e sites. For this reason, we believe that the phase transition from tetragonal to cubic is entropy driven. Li di usion is likely initiat ed at the 8a tetrahedral sites, requiring the neighboring 32g octahedral lithium to simultaneously di use into the unoccupied 16e tetrahedral site. However, lithium at 32g site mainly act as relay atoms and only contribute slightly to the direct filling of 16e site at high temperatures. Finally, we see that in both the tetragonal and cubic phase that lithium is locally distributed in clusters in only a few ways. These local clusters are mostly symmetrical - st of the di usion through the bottleneck. 92 C HAPTER 6 : Lithium self - diffusion in a model lithium garnet oxide Li 5 La 3 Ta 2 O 12 : a combined quasi - elastic neutron scattering and molecular dynamics study 6.1 Abstract I n this work lithium self - diffusion in the model lithium garnet oxide Li 5 La 3 Ta 2 O 12 was investigated by combining quasi - elastic neutron scattering experiments and molecular dynamics simulation. The Q - dependence of the quasi - elastic broadening measured experimentally and of the mean relaxation rate of lithium calculated by molecular dynamics simulations were both well described by the Singwi - Sjölander diffusion model. This model describes lithium nuclei that undergo diffusion via a jump - type mechanism consisting of a distribution of jumps with a mean square jump distance and residence time. The extracted mean jump length is consistent with a jump between tetrahedral (24 d ) to octahedral (48 g and 96 h ) sites within the Ia d symmetry structure and the residence time of Li at these sit es obeys an Arrhenius relation from ~ ps timescales at 1100 K to the ~ ns range at 400 K. This result supports a lithium diffusion mechanism in which jumps are more frequent (smaller residence time) and shorter at higher temperature. The self - diffusivities of Li from both experiment and calculation were in good agreement, but deviations from those previously measured using nuclear magnetic resonance and muon spin relaxation were observed and discussed. Analysis of stretching parameter describing the relaxat ion of lithium calculated by molecular dynamics indicated that Li motion is more cooperative at shorter length scales, below ~ 4.4 Å, which corresponds to the distance between octahedral sites across the tetrahedral site. 93 6.2 Results and d iscussion 6.2.1 Static structure factor T r eveals three Bra gg peaks ( Figure 30 (a) ) . For comparison, the static structure factor, , from the 700 K MD simulation is also presented in Figure 30 and is in agreement with that experimentally measured. The four sharp features are the 112, 022, 123, and 004 reflections from the Ia d structure 3, 37 . Experimental peaks contain instrumental broadening. 6.2.2 Overall intermediate scattering function and dynamic structu re factor T he neutron - weighted ISF calculated at 700 K is useful for understanding the components of the measured DSF that will be used to form a phenomenological model, and is presented in Figure 30 (b) . The ISF is shown (in later sections) to consist of vibrating La, Ta, and O with the only diffusing species being Li. Among the four species, only Li and La have Figure 30 : (a) The static structure factor obtained from the QENS experiment and MD simulations at 700 K. The results are vertically offset for clarity and lines serve as a guide to the eye. (b) Neutron - weighted (including both coherent and incoherent intermediate scattering functions calculated from MD at 700 K for several Q (a) (b) 94 significant incoherent scattering power ( Table 4 ) . I ( Q , t ) for vibrating species contributes to the plateau in Figure 30 (b) that becomes a delta function in the energy domain, S ( Q , E ), whilst the clear decay in I ( Q , t ) arising from Li becomes a Lorentzian - like function in S ( Q , E ) . Therefore, it is expected that the experimentally - measured S ( Q , E ) is well described by a linear function (arising from the background) B ( Q , E Q , E ), and a Lorentzian function (arising from diffusing Li) L ( Q , E ) , convoluted with the resolution function : Equation 43 where I 1 and I 2 are intensities. The centers of both delta and Lorentzian functions were fixed to be zero and analysis was performed for Q < 1.1 Å - 1 and at 1.55 and 1.65 Å - 1 to avoid coherent Bragg peaks ( Figure 30 (a) ), i.e. to limit the analysis to incoherent scattering and self - diffusion. The measured S( Q , E ) at 700 K for selected Q and at different temperatures for Q = 0.45 Å - 1 are shown in Figure 31 ( a ) and Figure 31 ( b ) , respectively. An example of the fit of the model to the data is shown Figure 31 ( c ) for 700 K at Q = 0. 35 Å - 1 . With increasing Q and temperature, the peaks broaden, with data at 300 and 400 K ( Figure 31 (b) ) close to that of the instrumental energy resolution (~ 3.0 µeV) due to limited Li mobility. A plot of the half - width - half - maximum (HWHM, ) of the Lorentzian function with Q 2 is shown in Figure 31 ( d ) . A linear dependence of with Q 2 ( DQ 2 , dotted line) is characteristic of continuous random walk 95 translational diffusion, also called Fickian diffusion, where D is the self - diffusion constant. The broadening of with Q 2 at low Q is linear and reaches a horizontal asymptotically at higher Q as is typical of a jump - diffusion process ( Figure 31 ( d ) ). The Q - dependent broadening of , describing the t ranslation of Li, was described well by a special version of the Singwi - Sjölander (SS) jump diffusion model 119, 138 in which diffusion occurs via jump diffusion where the time taken for the jump is much shorter than the residence time of the diffusing nucleus at a site, . This SS model is characterized by , where jumps are described by a spatial Figure 31 : (a) Dynamic structure factor S(Q, E) from QENS experiments at 700 K for different Q. (b) Dynamic structure factor S(Q,E) from QENS at Q = 0 .45 Å - 1, for different temperatures. (c) An example of the Wens fit with background, delta, and Lorentzian functions convoluted with the resolution function for 700 K and Q = 0.35 Å - 1, along with the residuals of the fit. (d) HWHM ( ) of the Lorentzian as a function of Q2 at different temperatures. Solid lines are the fits to the Singwi - Sj öl ander model. 96 probability (gamma) dist ribution with a mean square jump distance of . At low Q the model reduces to the Fickian model with . The residence time and mean jump length are 7.5 0.9 ps and 1.96 0.15 Å, respectively, at 700 K. The Q - dependent broadening of at 600 and 500 K is also shown in Figure 31 ( d ) , alongside a fit of the SS model. Figure 32 : (a) Intermediate scattering functions of species La, Ta, and O at 700 K from Md for selected Q. (b) fit of the plateau - Waller formula. (c) Comparison of the mean square displacement Uiso from MD with that obtained from neutro n powder diffraction. 3 97 6.2.3 Species - resolved dynamics: La , Ta, and O T he ISF for individual atomic species (La, Ta, and O) at 700 K calculated from MD simulations are shown in Figure 32 (a) for selected Q . A rapid exponential d ecrease in the ps range was followed by a stable non - zero plateau, which suggests that these species are vibrating instead of diffusing, as consistent with the Li 5 La 3 Ta 2 O 12 crystal structure ( Figure 32 ( a ) ) 3, 37 . For a harmonic oscillator, the structure factor has a Debye - Waller formula as , where is the mean square displacement (atomic displacement parameter) and the plateau values ( Figure 32 (a) ) were in good agreement with this. Similarly, the extracted are in good agreement with those obtained from neutron powder diffraction 3 as shown in Figure 32 (c) . 6.2.4 Li dynamics I n contrast to La, Ta, and O species, I ( Q , t ) of Li at 700 K from MD ( Figure 33 ( a ) ) exhibit decay that is faster at higher Q , indicating diffusion. Ignoring data for the first 10 ps that undergo ballistic collision and other short - time dynamics, the decay was described by a stretched exponential (Kohlrausch - Williams - Watts, KWW 139 ) function , where is the stretching parameter and the relaxation time . An equivalent value , mean relaxation rate, can be obtained from where is the Gamma function. Although the physical meaning of the stretching parameter is not well established, the deviation from a purely exponential decay, i.e. has been considered previously to arise from cooperativity 140 or coupling 141 of atomic motions. 98 For example, was called the coupling parameter in the coupling model by Ngai 141 . The behavior of as a function of Q ( Figure 33 (b) ) reveals that although at small Q , this value decays and converges to a plateau at higher Q for all investigated temperatures. We note that a similar behavior of was predicted by Mo de - Coupling Theory (MCT) 124 and observed in the dynamics of protons in supercooled water obtained from MD simulation 142 . This Q - dependence of could be understood by comparing the length scale associated with Q and the intrinsic cooperative length in the material. At small Q (longer distance), the interaction between any pair of two atoms is weak so their motion becomes uncooperative or not coupled ( ~ 1). At higher Q (shorter distance), the short - range interaction renders atomic motion cooperative or coupled ( < 1). Co nsidering coupling behavior of Li dynamics in Li 5 La 3 Ta 2 O 12, the pure exponential decay ( ) at small Q indicates that the motion of different Li atoms are uncooperative (not coupled) and Li motions couple ( < 1) at higher Q and plateau at Figure 33 : (a) intermediate scatteri ng functions of Li at 700 K from MD for selected Q. The solid lines are the fit of the KWW model to data after the first 10 ps. (b) Q - dependence of the KWW stretching parameter at different temperatures. (c) Schematic of the nearest - neighbor tetrahedral (p ink) and octahedral (blue) cages in lithium garnet oxides. Filled yellow circles are the 99 Q 2 2.0 Å - 2 , yielding a cooperative length of 4.4 Å that corresponds to the approximate distance between octahedral sites across the tetrahedral site ( Figure 33 (c) ). Figure 34 compares results from the experimental measurement using QENS and calculated MD simulations for Li motions in Li 5 La 3 Ta 2 O 12 at 700 K, with the Q - dependent broadening of which are both described well by the SS jump diffusion model. If we also plot from the QENS data ( Figure 31 ( d ) in Figure 34 we can see that diffusion is slightly faster from experiments than from MD simulation. We note that experimental data were analyzed in the energy domain using a Lorentzian function, which c orresponds to an exponential function in the time domain, while the analysis of MD results utilized a stretched exponential. We suppose this difference in data analysis might cause the mismatch, in addition to the intrinsic limitation of the empirical forc e fields. (a) Figure 34 : (a) HWHM ( ) experimentally measured using QENS and equivalent of Li diffusion derived from MD simulations in L5LT at 700 K. Solid lines are the fit of the SS model. (b) Self - diffusivity of Li in lithium garnet oxides of different compositions at different temperatures using various measurement probes. The activation energies for both MD and experimental sources are also provided. (b) 100 The residence time and jump length obtained from the SS model fits to both MD simulation and QENS experiments are shown in Figure 35 . These values are in good agreement, and the residence time of Li in the investigated temperature range is in the ps range, with slightly higher values obtained fr om MD simulation. Activation energies obtained from an Arrhenius fit are 0.34 0.02 and 0.28 0.03 eV for MD and QENS, respectively. These results are shown alongside those previously reported for the lithium garnet oxide Al - doped Li 7 La 3 Zr 2 O 12 , obtained from SLR - NMR experiments 54 , in which Li is reported to reside in the cage an order of magnitude longer than we find for Li 5 La 3 Ta 2 O 12 . We note that there are significant differences in the measurements used to obtain these results, with QENS measuring a density correlation and SLR - NMR measures the spin Hamiltonian correlation. The mean jump length obtained from QENS data is ~ 1.9 Å between 5 00 and 700 K, which is roughly the distance between 24 d (tetrahedral) and 48 g (octahedral) sites in Li 5 La 3 Ta 2 O 12 from previous results 3 . As the temperature increases to 1100 K (from MD), this jump length decreases to 1.5 Å, corresponding to the distance between 24 d (tetrahedral) and the closest off - center 96 h (octahedral) site, whilst at lower temperatures (400 K) the jump distance Figure 35 : (a) Mean residence time ( ) of Li diffusion in Li5La3Ta2O12 at sites obtained from MD simulation and QENS experiments. Solid lines are an Arrhenius fit. The residence time obtained from SLR - NMR experiments are also shown for Al - doped Li7La3Zr2O12 (cubic, C - L7LZ). (b) Mean jump length of Li diffusing in Li5La3Ta2O12 obtained from MD simulation and QENS experiments. Solid lines are guides to the eye. (c) Schematic of the relation between crystallog raphic sites for Li within tetrahedral (Td) and octahedral (Oh) cages in Li5La3Ta2O12 showing key distances from previous results. 101 increases to 2.5 Å, corresponding to the distance between 24 d (tetrahedral) and the f arthest off - center 96 h (octahedral) site. These results suggest a lithium diffusion mechanism in which jumps are more frequent (a smaller residence time) and shorter at higher temperatures, consistent with the delocalization of Li in Li 5 La 3 Ta 2 O 12 at higher temperatures previously measured 61 . The dynamics of Li in Li 5 La 3 Ta 2 O 12 that we measure and calculate lead to translational self - diffusion constants as shown in Figure 34 (b) . Since MD simulation results yield a similar mean jump length but slightly larger average residence time compared to those experimentally measured using QENS, the diffusivity determined from MD is lower than t hose obtained from the QENS measurement. The derived activation energies from the Arrhenius fit are 0.30 0.01 and 0.27 0.02 eV for MD and QENS, respectively, in good agreement with each other. The Li diffusivity obtained from other experimental probes for similar materials are also plotted in Figure 34 . These include PGSE - NMR measurements of Li 6.6 La 3 Zr 1.6 Ta 0.4 O 12 by Hayamizu et al 56 , of Al - doped Li 7 La 3 Zr 2 O 12 by Hayamizu et al 57 , and of Li 5.22 Al 0.26 La 3 Zr 1.5 W 0.5 O 12 by Wang et al 58 , as well as µSR measurements of Li 5 La 3 Nb 2 O 12 by Nozaki et al 59 and Li 6.6 Al 0.25 La 2.92 Zr 2 O 12 by Amores et al (we took the units discussed in the text, cm 2 /s, instead of what was labeled in the figure, m 2 /s) 60 . The fast Li diffusivity in Li 5.22 Al 0.26 La 3 Zr 1.5 W 0.5 O 12 was attributed to the measurement capturing only faster octahedral Li motions 58 . Diffusivity of Li in Li 6.6 La 3 Zr 1.6 Ta 0.4 O 12 and Al - doped Li 7 La 3 Zr 2 O 12 (Hayamizu et al 56 - 57 ) are slightly lower than we determine for Li in Li 5 La 3 Ta 2 O 12 , while those determined for µSR in Li 5 La 3 Nb 2 O 12 (Nozaki et al 59 ) and Li 6.6 Al 0.25 La 2.92 Zr 2 O 12 (Amores et al 60 ) are significantly lower. Compositional differences are likely to play a part in the se differences, and our future work will focus on studying Li diffusivity in Zr - doped Li 5 La 3 Ta 2 O 12 to understand these effects. Additionally, different probe methods capture different processes, with PGSE - NMR directly measuring 102 diffusivity at Q 0 and in th e s - ms time scale in which Li motions at the crystal surface and defects are sampled, and µSR measuring changes to the muon site local field and indirectly probing residence time. 6.3 Conclusions I n this work, we combine backscattering QENS experiments ( energy domain measurements) and MD simulation (time domain calculations) to study Li self - diffusion in Li 5 La 3 Ta 2 O 12 . Both QENS experiments and MD simulation suggests that Li motion follows a jump - diffusion type mechanism that is well described phenomenolog ically by the Singwi - Sjölander model in which Li move between crystallographic sites with a distribution of jump lengths . The Li motion becomes more cooperative at length scales shorter than ~ 4.4 Å, corresponding to the approximate distance between octahedral (48 g and 96 h ) sites across the tetrahedral (24 d ) site. The extracted residence times for Li at these sites obeys an A rrhenius relation to temperature. The extracted average jump lengths between 500 and 700 K are consistent with jumps between 24 d (tetrahedral) and 48 g (octahedral) sites . As the temperature increases to 1100 K the jump length decreases, corresponding to ju mps between 24 d (tetrahedral) sites and the closest off - center 96 h (octahedral) sites, whilst at lower temperatures (400 K) the jump length increases, corresponding to jumps between 24 d (tetrahedral) sites and the further off - center 96 h (octahedral) sites. This result supports a lithium diffusion mechanism in which jumps are more frequent (smaller residence time) with shorter at higher temperatures, consistent with the delocalization of Li at elevated temperatures. 103 CHAPTER 7 : Finite - size e ffects on the m o lecular d ynamics Simulation of f ast - ion c onductors: A case study of lithium garnet oxide Li 7 La 3 Zr 2 O 12 7.1 Abstract A useful tool to study ionic conduction mechanisms in fast - ion conductors is the molecular dynamics (MD) simulation performed on finite simulation cells with periodic boundary conditions. While there have been many studies of the effect of cell size on th e thermodynamics and kinetics of simple liquids, the finite - size effect in fast - ion conductors remains elusive. This work presents a case study to investigate the finite - size effect on the phase transition, self - diffusivity, ionic conductivity, Haven rati o, thermodynamic factor, and Fickian diffusivity using lithium garnet oxide Li 7 La 3 Zr 2 O 12 as a model material. It was found that cell sizes influence extracted thermodynamic and kinetic characteristics in different ways with magnitude ranging from weak to strong. For Li 7 La 3 Zr 2 O 12 , reliable properties can be obtained with a 3×3×3 (5184 atoms) cell. 7.1.1 Background M olecular dynamics (MD) simulation s for crystalline materials are generally performed within a finite simulation cell with periodic boundary con ditions. In many MD simulation studies, it was observed that structural and dynamical properties obtained from MD results varied with the size of the simulation box. For example, the self - diffusivity of liquid atoms was often found to increase with incre asing system size. 143 - 146 In addition, small subsystem s are usually carved out from the cell to simulate the grand canonical ensemble and extract thermodynamic properties. These thermodynamic properties were found to depend strongly on the size of subsystems. 126, 130 However, most investigations into these finite - size effects have been centered on simulating water or Lenard - Jones fluids while leaving fast - io n conductors largely unexplored. 104 7.2 Results and discussion 7.2.1 Phase characterization A n important aspect in determining the validly of our simulations is to ensure that they accurately reproduce lattice parameters and phase transition characteristics r epresented experimentally. In the literature, a wide variety of transition temperatures ranging from 373 to 923 K was reported. 4, 30 - 31, 43, 98 Literature results and our previous work suggest that such variation could be attributed to intentional or inadvertent impurity (e.g. Al)/dopant incorporation, H 2 O or CO 2 exposure, all of which reduces the transition temperature. 1, 30, 98 Thus we chose experimental data with the highest transition temperatures, i.e. works of Larraz et al 43 and Matsuda et al 98 , as the benchmark. Figure 36 (a) shows the latt ice parameters (normalized to 1×1×1 cell) for different cell sizes resulting from the sequential cooling from 1400 K, along with experimental measurements. Overall, our simulation results match closely to experimental ones in terms of change of lattice pa rameters and phase transition temperature (~900 K), with the exception of 2×2×2 cells indicating a transition temperature of 1000 K. The only major deviation in average lattice parameter shown is for the 1 ×1×1 simulation at low temperatures, where it cons istently has larger lattice parameters compared to the larger simulations . 105 (a) (b) (c) Figure 36 : Finite - size effects on the phase transition. (a) Average lattice parameters as a function of temperature for 1×1×1 4×4×4 simulations. Solid lines are guide to the eye. Experimental data from Matsuda et al 98 and Larraz et al 43 are shown for comparison. (b) Fluctuation of lattice paramet ers as a function of number of atoms for selected temperatures. Solid lines are the linear fit (slopes close to - 0.5). (c) Time evolution of lattice parameters for different cell sizes at 900 K. The standard deviation of lattice parameters decreases wit h the increase of cell sizes in Figure 36 (a). These values are plotted in Figure 36 (b) as a function of number of particles to further investigate the size effects. For simplicity, we only plotted data for 800, 1100, and 1400 K. For the ideal gas with number of particles and cubic length, the fluctuation of in a NPT ensemble takes the form of . The linear fit of the log - log plot in Figure 36 (b) yields slopes varying from - 0.49 to - 0.50, consistent with the dependence of ideal gas. Another way to study size effects is to examine the fluctuation of lattice parameters around the phase transition when cells are volatile. Lattice param eters (as orthogonal a, b, c) are shown as a function of time at 900 K in Figure 36 (c). Larger 3×3×3 and 4×4×4 simulations both show very small variat ion in lattice parameter and maintain a mostly cubic phase throughout the simulation. The 2×2×2 simulation shows a more complicated situation where a slightly tetragonal unit cell undergoes reorientation of its lattice vectors as the simulation completes. We treated 2×2×2 simulation as a tetragonal phase and extracted the average lattice parameters from 106 25 - 60 ps. The 1×1×1 simulation shows larger fluctuations with more frequent tetragonal reorientation, making an accurate description of phase difficult. To simplify the extraction of lattice parameters, we characterized it as being cubic. After determining that 3×3×3 and 4×4×4 cells were large enough for the phase characterization, we further explored the hysteresis of phase transition in L7LZ by sequentially cooling from 1000 K and heating from 800 K, respectively, in 25 K steps. Figure 37 plots the lattice parameters in a heating or cooling cycle. For 3×3×3 cells, we see that starting at 800 K and heating, the tetragonal to cubic transition occurs at aro und 950 K, but during cooling we did not observe a tetragonal phase until below 875 K. For 4×4×4 cells, the hysteresis region was slightly reduced to within around 25 K. Similar hysteresis behavior was also observed in our previous ionic conductivity mea surements. We think 4×4×4 cells capture the hysteresis more accurately but 3×3×3 cells offer a good balance of accuracy and efficiency. Figure 3 7 : Change in lattice parameter for 3x3x3 and 4x4x4 under heating and cooling cycles 107 7.2.2 Self - diffusivity, ionic conductivity, and Haven ratio Self - diffusivity rep resents the ability of an individual atom to move. The self - diffusivity values of Li ions for the four cell sizes are shown in Figure 38 . Our calculated values are sligh tly higher than those from Jalem et al 47 employing first - principles molecul ar dynamics (FPMD) simulation with 1×1×1 cell sizes. The difference could be caused by the difference in lattice parameters employed. Experimental measurement of Li self - diffusivity by the spin - lattice relaxation Nuclear Magnetic Resonance (NMR) studies yielded diffusivities on the order of 10 - 14 cm 2 s - 1 at ~300 K. 31 Our simulations were unable to calculate diffus ivity at such a low temperatures due to the low mobility of Li, but extrapolating from 800 - 700 K for the 3×3×3 and 4×4×4 simulations yields a diffusivity on the order of 10 - 17 cm 2 s - 1 . However, it is worth noting that two relaxation times that differ by ar ound 3 orders of magnitude were reported in the NMR study and the faster relaxation time was employed to calculate the diffusivity. If the slower Figure 38 : Self - diffusivity as a function of temperature for four different simulation size s. Literature data from the FPMD simulation (1×1×1) by Jalem et al are shown for comparison. The inset shows details at high temperatures. Solid lines are guide to the eye. 108 relaxation time was used instead, the result would be closer to the value obtained in the present computation al work. Previous MD simulations have shown that self - diffusivity of simple liquids such as water 143 and silica 144 increases with increasing cell size in the form of , where is the diffusivity at cell size , the diffusivity at the thermodynamic limit, the Boltzmann constant, the temperature, and the shear viscosity. While the present study indicates that 1×1×1 structures maintain a consistently higher diffusivity than that in larger simulation cell s, overall we observe no size - determining trend for 2×2×2, 3×3×3, and 4×4×4 simulations in Figure 38 . Such a lack of clear size effects could be due to dif ferent interaction potentials in L7LZ compared with those in simple liquids. With respect to the appropriate simulation size, the 2×2×2 simulations allow for a reasonable calculation of diffusivity at higher temperatures, but as we approach the transition temperature the increased stability of the larger 3×3×3 and 4×4×4 cells allow for a more reasonable diffusivity calculation. Next , we will study the size effect on the ionic conductivity, one of the most important properties of fast - ion conductors. The io L7LZ has been widely studied, and is commonly used to draw comparisons between experimental and simulation properties. 29, 34, 38 - 39 For the reason discussed earlier, we selected the experimental conductivity results with high transition temperatures (~900 K), e.g. from Wang et al 1 and Matsui et al 4 . For the computational results, we first applied the Nernst - Einstein relation for ideal solutions, i.e. Equation 8 , to convert self - diffusivity to conductivity, as commonly done in the literature, e.g. in the classical MD simulation by Adams et al 50 and FPMD simulation by Jalem et al 47 and Miara et al. 5 The computed conductivities for 3×3×3 cells are s imilar to those from Miara et al 5 (1×1×1 cell) in Figure 39 . For nonideal solutions such as L7LZ , we plot the 109 ionic conductivity calculated from Equation 19 in Figure 39 for all cell sizes. Similar to what was observed for self - diffusivity, ionic conductivity shows small variations fr om 2×2×2 to 4×4×4 cells. The difference between calculated and experimental values in diffusivity and conductivity could be due to impurity incorporation of the experimental samples, e.g. Al. Such impu rity effect will be reported in future publication s . We are also in the process of exploring other means of calculating the conductivity with respect to the methods outlined in Table 4 . It appears out model may be under - binding for lithium despite the good agreement with structural parameters. As discussed in the Chapter 3 , the Haven ratio characterizes the ratio of self - diffusion vs collective diffusion (ionic conductivity). Its value is commonly assumed to be 1, i.e. no correla tion of motion of different atoms. This allows the use of self - diffusivity and Nernst - Figure 39 : Finite - size effects on the ionic conductivity for 3×3×3 and 4×4×4 cells. Values calculated from both Equation 8 Equation 19 are presented. Literature data from experimental measurements by Wang et al 1 and Matsui et al 4 , along with FPMD simulations by Miara et al 5 110 However, for complex materials such as fast - ion conductors, the Have n ratio is expected to deviate from 1. Finite - size effects on the Haven ratio are shown in Figure 40 . At high temperatures, Have n ratios are around 0.3 and size dependence is weak. At lower temperatures, the ratios dropped to a few percent signifying strong collective motion. 7.2.3 Thermodynamic factor and Fickian diffusivity I n fast - ion conductors, the thermodynamic factor reveal s the effective fraction of mobile charge carriers, as discussed in the Appendix. There are three methods in obtaining the thermodynamic factor from the trajectory of MD simulation: fluctuation, KBI, and static structure factor. The fluctuation method involves dividing a simulation cell into subsystems that are large enough compared with the microscopic volume but small enough compared with the cell volume ( to achieve grand canonical ensemble). 126, 147 We applied this method in our previous study of Li 7 La 3 Zr 2 O 12 by averaging values corresponding to of 1/8 and 8/125 for 2 × 2 × 2 cells. 148 How ever, the guideline to choose the right size of Figure 40 : Finite - size effects on the Haven ratio for four simulation cell sizes. Solid lines are guide to the eye. 111 subsystem remains unclear. The KBI method involves evaluating the Kirkwood - Buff integral of the Li Li pair distribution function 130 . For simple liquids in which the pair distribution function decays to the value of one quickly, this method allows reliable estimation of . For complex materials such as fast - ion conductors, such KB inte grals are more susceptible to finite size truncation errors. An alternative is the static structure factor method which is the KBI method in the Fourier space. This is the method employed in the present study. - 1 from different cell sizes at different temperatures are shown in Figure 41 . Overall, with increasing temperature we see a - 1 suggesting that L i diffusion is becoming less coordinated at higher temperatures. A kink at 900 K was identified for all cell sizes, corresponding to the tetragonal to cubic phase transition. Since the thermodynamic factor is related to the number fluctuation reflecting interatomic interaction, the two regions above and below 900 K imply different Li Li - 1 values are Figure 41 : Fin ite - - 1 at different temperatures. Solid lines are guide to the eye. 112 decreasing with the increasing cell sizes and they start to converge above 3×3×3 cells. The - 1 in both 3×3×3 and 4×4×4 cells are on the order of a few permille which suggest that only a small fraction of Li ions contribute to the conduction. A value of ~6% was reported in the literature based on the muon - spin relaxat ion and conductivity measurement for L7LZ at 300 K, although the reported conductivity is similar to our calculated value at around 800 K. 59 Such difference in effective concentration is likely caused by the impurity incorporation discussed earlier for the ionic conductivity. Finite - size effects on the Fickian diffusivity are shown in Figure 42 . According to Figure 42 , the Fickian diffusivity is related to the ionic conductivity ( Equation 19 ) through the thermodynamic factor ( Figure 41 ). While the ionic conductivity has weak dependence on the cell size, the Fickian diffusivity increases with increasing cell sizes due to the opposite trend in t hermodynamic factors. Figure 42 : Finite - size effects on the Fickian diffusivity for four simulation sizes . 113 7.3 C onclusions I n this work, we performed the molecular dynamics simulation and investigated the effect of simulation cell sizes on the phase transition, self - diffusivity, ionic conductivity, Haven ratio, thermodynamic factor, and Fi ckian diffusivity of a model fast - ion conductor: Li 7 La 3 Zr 2 O 12 . In terms of phase transition, even the smallest cell size (1×1×1 with 192 atoms) was able to capture the tetragonal to cubic phase transition. Increasing the cell size reduces the fluctuation of lattice parameters in the form of and increases the accuracy in estimating the phase transition hysteresis. No obvious size effect was detected for the self - diffusivity of Li 7 La 3 Zr 2 O 12 , as opposed to the observation of increasing diffusivit y with increasing cell size for simple liquids. The ionic conductivity and Haven ratio also shows similar weak size dependence as the self - diffusivity. In terms of thermodynamics, the inverse thermodynamic factor decreases with the increasing cell sizes, with values on the order of a few permille suggesting only a small fraction of all lithium - ion s are responsible for the electrical conduction. On the contrary, Fickian diffusivity increases with increasing cell sizes since it is inversely related to the inverse thermodynamic factor. Overall, a cell size of 3×3×3 (5184 atoms) offers the best balance between accuracy and efficiency for Li 7 La 3 Zr 2 O 12 . The methods and consideration employed in this study can be extended to other fast - ion conductors. 114 CHAPTER 8 : Modeling Li 7 - xLa 3 Zr 2 - x Ta x O 12 using core only model 8.1 Abstract C ritically, one of the most important areas of investigation regarding the lithium garnet series is the effect of li thium concentration coupled with the tantalum and zirconium substitution. Based upon the results from the finite size effects on L7LZ, 3 x 3 x 3 simulation cells were generated for Li 7 - x La 3 Zr 2 - x Ta x O 12 for x = 0, 0.25, 0.5, 0.75, 1, 1.5, & 2. Within this framework we were able to show that the phase transf ormation temperature is a function of composition , where the critical lithium concentration of 6.5 per formula unit (p.f.u) is observed with agreement to reports by Thompson et al 42 for aluminum substitution. Our simulations predict a maximum c onductivity for concentrations of Li in the range of 6 .5 to 6. 7 5 p.f.u. in agreement with previous reports . 35 The inverse thermodynamic factor for Li Li interactions ( LiLi - 1 ) is maximized for a lithium concentration of 6 p.f.u. implying this composition allows for the highest degree of mobile species. for nonidentical pairs , with respect to Li and the frame work atoms La, Ta, Zr, and O are also investigated. We show that the framework atoms may play a crucial role in the redistribution of lithium across the Li sublattice maximizing the entropy of the system , and encouraging higher self - diffusivity values when the con centration of Li is around 6.5 p.f.u. 115 8.2 Results and d iscussion 8 . 2.1 Phase characterization of 3x3x3 LLZT cells P revious work show ed that 3 x 3 x 3 simulations a re sufficiently large and efficient enough to model the phase transformation in L7LZ . Continuing upon that work, potential set verification was performed by comparing the simulated structural properties of the LLZT series to X - ray and neutron diffraction v alues in the literature. Lithium garnets of composition have consistently shown the cubic single - phase behavior experimentally and the anisotropic NPT simulations performed on these compositions concur with reported values and are plotted in Figure 43 (a) . Similarly, when comparing the volume of unit cells within this compositional range to X - law is observed, in accordance to other simulation studies. There are however two noticeable devi ations. First the predicted volume of L5LT (x=2) in Figure 43 (b) is about 1.4% less than the values reported from experiments. Second, both X - ray and n eutron studies have shown that there L7LZ which is not captured by most simulation Figure 43 : Anisotropic NPT simulations on 3x3x3 simulation cells. (A) Unit cell lattice parameters for LLZT series. Solid lines are guides to the eye. (B) Unit cell volume at 300 K, 600 K, 900 K, and 1200 K for LLZT series. Experimental data from Logeat et al, Mukh opadhyay et al, and Wang et al. are shown for comparison. (C) Time dependent lattice parameter for L - 6.75LZT at 300 K, 400 K, 500 K, 600 K. (a) (b ) (c ) 116 studies. It is likely that the loss of volume conservation is related to impurity incorporation from moisture exposure causi ng lattice distortions and the creation of multiple garnet compositions. For compositions, the simulations suggest that at low temperatures the tetragonal phase is preferable. Figure 43 (c) plots orthogonal a, b, and c lattice vectors as a function of time. At 500 K and 600 K the a, b, and c lattice remain relatively constant while fluctuating within 0.1 Å of the mean. As the temperature is lowere d to 400 K, a separation begins to form between the lattice vectors as the cubic phase transforms to a slightly tetragonal phase by the end of the simulation. Finally, when the simulation cell is cooled to 300 K we see a consistent tetragonal phase indicat ing that this composition of garnet may be unstable. There is some evidence that lithium will self - separate into tetragonal L7LZ and cubic L 6.5 - y LZT. 149 In compositions where the critical Li concentration has been exceeded leading to evidence that a single - phase solution is not formed but instead a mixture of L 6.5 LZT and L7LZ . Due to the relatively small simulation cell used in this study it is not possible to model the possib le formation of multiple phases, so further discussion on structural properties of L 6.75 LZT will assume a single tetragonal phase is formed. From this assumption if we consider L7LZ to be the starting composition, Ta substation on the 16a sites works in a similar manner to other dopants. When Zr is substituted for Ta the cubic to tetragonal phase transition temperature is depressed. For compositions above the critical Li concentration ( ) the tetragonal phase is completely suppressed. In the case of L 6 .75 LZT, simulations suggest that the cubic to tetragonal transformation it is depressed to ~500 K versus ~900 K for uncontaminated L7LZ . Excess entropy calculations discussed in a subsequent section of this paper suggest a more ordered sub - lattice configur ation reminiscent of tetragonal L7LZ is present in L 6.5 LZT making a definitive answer on the stability of this composition difficult to determine. 117 8 .2 .2 Ionic conductivity, T he collective movement of Li ions denoted as the ionic conductivity has been proven to exhibit unique phenomena related to the stoichiometry of Ta and Zr at room temperature. The expected trend of having the lowest conductivity at the lowest Li concentrat ion ( L5LT , x = 2) and highest conductivity at the highest Li concentration ( L7LZ , x = 0) is disrupted with the actual maxima disputed to be in the range of depending on the quality of the sample. The dramatic decrease in conductivity in L7LZ results from the ordering of Li on the occupied 8a, 16f, and 32g sites and the unoccupied 16e sites. Figure 44 (a) plots Li conductivity from 1400 K to 500 K with an extrapolation down to 300 K. Except for L7LZ each composition has a near linear profile in accordance with Arrhenius behavior. Upon closer inspection, it is evident that around 900 K, each composition exhibits a slight kink in their slope pr ofile. This kink has can be observed with most garnet compositions experimentally at around 500 K and is potentially a sign of a glassy transition of the Li sub - lattice. (a) (b ) (c ) Figure 44 : (a) Arrhenius plot of conductivity for LLZT series. Squares are calculated conductivity values; lines are best fit from 1000 K 500 K then extrapolated down to 300 K. Experimental total conductivty for LLT shown by Thangadurai. Inset shows high temperature conductivity with lines as guides. (b) Iso - temperature plot of conductivity as a function of lithium concentration. 1200 K, 900 K, 600 K are direct measurements, 300 K is extrapolated from best fit line in (a). (c) Activation energy as a function of lithium concentration for LLZT series Co mpared to reported activation energies for various garnet compositions with standard deviation bars. 118 Figure 44 (b) plots the conductivity of Li as a function of Li concentration. Due to poor kinetics at 300 and 400 K, the values shown at 300 K are that of the extrapolated line from Figure 44 (a). The highest room temperature conductivity observed in the simulation is that of Li 6.75 LZT, but because this data set is extrapolated from its cubic phase the effects of having a tetragonal phase are not captured in this analysis. Based upon sub - lattice structural effects described later in this work it is likely that L 6.5 LZT is the most conductive phase at room temperature. The experimental values by Wang 35 indicate good agreement with the simulations where . Their reported maximum being located at Li 6.75 LZT is likely from Al contamination and stabilization derived from the synthesis and sintering process . I t is likely the a ctual concentration of Li in the bulk is less than the theoretical stoichiometry. From the best fit lines in Figure 44 (a), the activation energy (E A ) for each composition was calculated and plotted in Figure 44 ( c). L5LT to L 6.25 LZT have a mainly flat profile with a slight peak and valley at L 5.5 LZT and L 6 LZT respectively ranging from 0.34 to 0.31 eV. When the concentration of Li increases from 6.25 to 6.75 there is a decrease in activation energy equal to about 0.1 eV to a value of approximately 0.22 eV before shooting up dramatically for the tet ragonal L7LZ . When compared to the activation energy for all Li garnet type compositions in the review by Thangadurai 34 , there is a strong agreement in the general trend, but a difference of about 0.15 eV for all measurements 33 . The d iscrepancy between observed activation energies and those modeled is likely resulting from reports of total resistivity which includes a merging of the bulk and grain boundary signals. Typically, these will increase the activation energy and they a feature that is not modeled in my system. 119 8.2 .3 Self - diffusivity While the conductivity discussed above concerns itself with the collective motion of Li atoms the self - diffusivity allows for discussion about the ability of individual Li atoms to jump between neighboring sites. Ideally, the two properties should be highl y correlated. In the case of an ideal gas where the Nernst - Einstein relation is valid, the motion of each species is independent of its neighbors. While this assumption has been used quite liberally previously, lithium garnets have a relatively small conce ntration of independent mobile species and instead primarily move through collective motions making calculations of conductivity through the Nernst - Einstein relation incomplete without a corrective term. Self - diffusivity and conductivity exhibit nearly id entical trends regarding the LLZT series. Evident from Figure 45 (a) and Figure 45 (b) as the concentration of lithium is increased up to 6.75 there is a corresponding increase in self - diffusivity terminating with a sharp decrease at L7LZ . A deviation is observed from this trend the x = 1.5 composition L 5.5 LZT, which has a lower self - diffusivity than that of the x = 2 composition L5LT . Except for L7LZ , which is much (a) Figure 45 : (a) Arrhenius plot of self - diffusivity for LLZT series with extrapolated best fit lines to 300 K. Inset shows high temperature 1 4 00 K to 1000K with lines acting as guides to the eye. (b) Isothermal plot of self - diffusivity as a function of lithium concentration. 300 K data taken from extrapolated lines of (a). ( b ) 120 lower due to the tetragonal phase transformation, the simulations predict self - diffusivity values four to six orders of magnitude lower the n what is seen with solutions of LiPF 6 in ethylene carbonate (EC) and diethyl carbonate (DEC). Despite the dramatic difference in self - diffusivity, the conductivities of EC/DEC solutions and that of LLZT are only one to two orders of magnitude apart. Makin g it evident that low a self - diffusivity for structured vacancy materials such as LLZT is not a major hindrance to the overall conductivity of ions. 8.2 .4 Haven ratio, thermodynamic factor T he following discussion will explore two measur es of collective and independent diffusion in the garnet series . Firstly, the Haven ratio ( H R ) expresses the ratio of self - diffusion ( vs collective diffusion ( of a species. Figure 46 (a) shows that the Haven ratio for LLZT series is in general increasing with temperature and lithium concentration. In the case of 300 K the data is quite noisy due to the poor kinetics at low temperatures. N evertheless , the low temperature trend of maximi zing the ratio of self to collective diffusion at L 6.5 LZT is observed up to 600 K where a maximum of about 0.3 is observed for all lithium concentrations . As the temperature is further increased, we see that for all compositions the Haven ratio begins to l evel (a) ( b ) Figure 46 : (a) Haven ratio for Li diffusion in the LLZT series at 300 K, 600 K, 900 K, and 1200 K and (b) for Li Li pairs at 300 K, 600 K, 900 K, and 1400 K . 121 out . E ffectively showing that increasing the l ithium concentration to 6.5 allows for more independent diffusion events at lower temperatures , essentially maximizing the ratio of independent diffusion to collective diffusion at a lower temperature comp ared to other compositions. It is possible that this composition is on the boundary of stability for the local cluster arrangements. In Chapter 4 we showed that the dependence of concentration on what types of cluster formed for L5LT and L7LZ respectively . Assuming that increasing the Li concentration is going to fo rce more Li - Li type interactions, it is possible that th is composition will undergo more local rearrangements than th e others, leading to a larger diffusivity. The second method for probing collective and independent motion , is the inverse thermodynamic factor ( ) . Here in contrast to the previous reports in Chapter 4 , it is evaluated from the static structure factor S ab ( Q , 0) as we approach the long - The respectively . If the correlated pairs are all the mobile l ithium type, when this condition is met the effective fraction of charge carriers is determined. Whereas in the case , species and the lattice are evaluated. For the sake of discuss ion, we postulate that only interactions consisting of the highest order of magnitude will be considered influential to lithium diffusion. Contrasting the single maximum in the conductivity, self - calculations, evaluation of presents a bimodal distribution. The two peaks in Figure 46 (b) located at L 6 LZT and L 6.5 LZT respectively are orders of magnitude greater than values for the end members L5LT and L7LZ having effective carrier concentrations of only a few permille. In 122 the exemplary case of an ideal gas , meaning all l ithium contribute to the conduc tivity, while no garnet composition approaches this limit, at all investigated temperatures the fraction of charge carriers approach ~11% in L 6 LZT. Increasing l ithium concentration further decreases the carrier concentration at 300 K in both L 6.25 LZT and L 6.5 LZT to ~6% & ~7% respectively. For the last intermediary composition L 6.75 LZT the fraction of carriers is further decreased to only ~1.5% at room temperature. Inferring from the data presented Figure 45 and Figure 46 , the increase in self - diffusivity associated with high l ithium concentration compositions plays a greater role in the resulting conductivity than the effective concentration of charge carrier s. I nterestingly, as the temperature increases for L 6.25 LZT and L 6.5 LZT there is a slight decrease in the fraction of effective charge carriers. Changes in have previously been related to changes in how lithium is dispersed about tetrahedral an d octahedral sites and the resulting lithium sub - lattice. L 6.75 LZT for example undergoes a change in similar to that of L7LZ . At temperatures sufficient to undergo the tetragonal to cubic phase transformation there is a noticeable increase in t he number of effective charge carriers. In a subsequent section, we discuss a descriptive reasoning for these results as it relates to the excess entropy of the l ithium sub - lattice. Now that Li Li interactions have been modeled the effects of the surround ing framework atoms need to be explored. Figure 47 plots the thermodynamic factor of framework atoms with mobile lithium species. For all end member compositions is relatively close to zero. This is indicative of a homogenous bulk s tructure where at the potential field around all pairs is nearly equivalent. As heterogeneity is introduced through Ta/Zr substitution what we call the effective potential field becomes unbalanced and a preferential pair coordination environment is forme d. With respect to Li - O and Li - La pairs, it can be assumed that these 123 framework ions have no effect on lithium conduction because is an order of magnitude lower than of Ta and Zr and close to zero. For Li - Ta and Li - Zr pairs however we can see that these ions act as a mirror to one another at each composition, Zr with its 4 + valance state correlates positively with lithium - ion s and indicat ing a preference between th e pairs . Ta with its 5 + state is negatively correlated with lithium positions over time suggesting t he imbalance in charge between these two atoms on 16a sites may act as a force for driving for the self - diffusion of l ithium . As more disorder is introduced the lithium may rearrange them sel ves in order to obtain preferential interactions. The Figure 47 : for (a) Li - La pairs, (b) Li - Ta pairs, (c) Li - O pairs, and (d) Li - Zr pairs . (a) ( b ) (c ) ( d ) 124 magnitude of at each composition likely correlates with the degree of heterogeneity on 16a sites. L 6 LZT with its equally numbered and randomly dispersed Ta and Zr species could allow for the most favorable distribution of Ta/Zr ensuring that the supporting f ramework of the lithium sublattice is sufficiently disordered. 8.2 . 5 Li site occupancy and regular solution model E xperimental neutron diffraction studies have shown a shifting in lithium site preference , predicting that with increased Li concentration, T d sites are depopulated as Oh sites become favored. 37 This has been reasoned to be a result of minimizing Li Li interactions , and a change in relat ive site energy. Because of the wide range of possible lithium distributions, and the highly cooperative diffusion mechanism involving the motion of nearest and next nearest neighbor dynamics it is infeasible to model each configuration and predict change s in site energy and how a given initial state will result in a final state using techniques such as nudge elastic band . In an NVE simulation however , it is possible to probe a wid e range of configurational e nergy and use the statistics of site occupancy an d average potential energy to form a regular solution model. First looking at the occupancy for both tetrahedral and octahedral sites for the garnet serie s. At the end member composition of L7LZ two phase behavior exists as 33% of Td sites are occupied at 300 K and 600 K accompanied by 100% occupancy of Oh sites. Once the simulation is heated and the phase transformation occurs , the total Td occupancy is increased to ~40% while Oh occupancy is decreased to ~96%. This agrees with our previous work on L7LZ using the core shell model and matches the trend of migrating some of the fully occupied 8a lithium to the unpopulated 16e sites giving a uniform distribution across all Td sites. 125 For other Li compositions, it is evident that increases in temperature cause a depopulation of 24d sites while increasing Oh occupancy. The regular solution model is constructed by assuming the 24d and 48g sites are separate sub - lattices with the interacti on energy between them defined by a mean - field term where the entropy of each sub - lattice is treated as if it is an ideal solution. A summary of terms is defined in Table 7 . Where d is the tetrahedral site occupancy, g the octahedral site occupancy, k B the Boltzmann constant, T the absolute temperature, the interaction energy, and 24d and 48g being the site energy of 24 d and 48 g sites respectively . Table 7 : Parameters for the regular solution model Number of 24 d sites, N ; constraint: 3 d +6 g = 5 + x x=[Zr] Site /Interaction Site Energy Entropy Interaction Energy 24 d 24d k B TN [d ln d + (1 - d) ln (1 - d)] 0 48 g 48g k B T2N [g ln g + (1 - g) ln (1 - g)] 0 24 d - 48 g (a) ( b ) Figure 48 : (a) tetrahedral and (b) octahedral occupancy as a function of Li concentration and temperature. 126 The term - La, Li - O, and Li - Ta/Zr pairs as follows, Equation 44 Combining all terms in Table 7 we are left with an expression for the free energy: Equation 45 The resulting expression is similar to that of a regular solution or Frumkin adsorption isotherem on a single lattice. For each simulation, the time averaged potential energy from NVE trials is taken to be F, while the corrective input term F intrinsic accounts for the configurational energy contribution from the framework atoms at 0 K taken from GULP simulations . The terms , and were then refined using the Td and Oh occupancy values from 0 K to 1400 K. The results of the refinement procedure are displayed in Table 8 . From the regular solution model, we can now see that as we transition from low to high Li compositions the relative site energy difference between 24d and 48g sites changes from slightly positive to strongly negative values implying that Td occupancy is fa vored at low Li concentration while Oh occupancy is favored at high concentrations. With respect to , negative numbers represent the degree to which Li is willing to coordinate with neighboring Li species. As the concentration of Li goe s up there is a consistent trend of decreasing interaction potential allowing for a greater nearest neighbor occupancy. The magnitude of for compositions is quite small compared to higher Li compositions suggesting that both lithium s ites are at similar relative energies and the driving 127 force for dispersed Li species stems from interaction parameters that suggest week coupling of Li Li pairs. For compositions with , octahedral sites become heavily preferred with values of implying that highly coordinated and occupied nearest neighbor clusters are formed. Overall this regular solution model was able to adequately fit our observed occupancy trends and provides insight on the degree of site preference for the LLZT series. Table 8 : Refined regular solution parameters [Li] (p.f.u) 48g - 24 d (eV unit - cell - 1 ) (eV unit - cell - 1 ) R 2 5 (x=2) 5.984 - 1.173 0.9973 5.5 (x=1.5) 8.307 - 2.787 0.9942 6 (x=1) - 4.688 - 6.415 0.9825 6.25 (x=0.75) - 20.182 - 13.12 0.9902 6.5 (x=0.5) - 63.158 - 29.97 0.9888 6.75 (x=0.25) - 180.53 - 72.18 0.9857 8.2. 6 Li sub - lattice and configurational energy T hus far we have examined Hr , , and Td/Oh occupancy values in pursuit of finding a descriptor for predicting what properties mark a highly conductive Li garnet composition. A final method to be explored is through the calculation of the excess entropy of the Li sub - lattice. Shannon 132 described the amount of information a system contains is proportional to the entropy of the system. From Equation 41 , where S is the entropy, k B is the Boltzmann constant, and p is the probability that a location within the unit cell is occupied is derived. To perform this calculation, we record the position of each atom for every frame of a 128 simulation and superimpose its position upon a three - dim voxels with uniform dimensions equal to approximately 0.1 Å. For every frame in which the pixel is occupied a value of 1 is assigned it. We then sum the pixels and divide by the total number of frames to find the probability of each pixel being occupied. In the case of an ideal gas each pixel will have a uniform probability for which the ideal entropy (S I ) of the system is maximized. The Li sub - lattice for the garnet series however has a finite number of positions in which Li can occupy, we call this the configurational entropy (S C ) for Li. The excess entropy for the system is then the difference of S I - S C . Because each of our unit cells is of varying length we then normalize the excess entropy by dividing by S I to ensure a uniform pixel density. Figure 49 plots the normalized excess entropy as both a function of composition (a) and a function of temperature (b) for the garnet series. We can see from the Figure 49 (a) that at room temperature the excess entropy is maximized for the [Li] = 6.5 composition. Upon heating this peak shifts to higher lithium concentration indicating that if the same amount of energy was added to the Li sublattice, more disorder is generated in the [ Li ] = 6.75 case vs the others . At 7 00 K, every composition except for L7LZ is maximally disordered for a given temperature , this change does not occur L7LZ unt o 900 K consistent with the tetragonal to cubic phase transformation. 129 The maximizing of the normalized excess entropy , we believe is proportional to the collective diffusion of lithium. When comparing the density maps of these garnet compositions across various temperatures we observe that as the concentration of lithium increases so too d oes the interconnectedness of our density maps. This type of response is similar to an idea proposed in the Adam - Gibbs entropy model of thermally activated transitions. 150 In the Adams - Gibbs theory , a change in activation energy is proportional to the inverse of the configurational entropy . Therefore, a composition with the largest normalized excess entropy (most disorder) will have a lower activation energy compared to the other compositions . If the entropy of the system can be defined rigorously it may be possible to modify our extrapolations of conductivity to lower temperatures. One of the major problems still being add ressed , is calculating low temperature properties where the cost of simulation time inhibits our ability to gather a statistical representation of diffusion events . The calculation of excess entropy is less computationally expensive than other methods and still valid if only vibrational or jump - back type motions are observed when an atom moves to one site and quickly returns to Figure 49 : Normalized excess entropy for the LLZT series as a function of composition and temperature. (a) ( b ) 130 the original . Such motions are difficult to capture in the mean squared displacement and auto correlation calculations, though they still represent low temperature events. By modeling how entropy changes with temperature it may be possible to correct for any non - Arrhenius behavior. As reported previously , we believe the transformation to the cubic phase at hig h temperatures appears to be driven by the change in entropy of the system. We can see a discontinuity in excess entropy between 800 K and 900 K , indicative of a first order phase transformation with respect to the dimensionless heat capacity of the system . The transformation coincides with the increase in conductivity reported previously in section 8.2.2 suggesting that increasing lithium disorder is the major contributor to the existence of a highly conductive phase. To visualize the changes in entropy we again can look to the density maps as a function of temperature for the various compositions. There are two major trends to be addressed, first the change in density for a single composition at different temperatures , and se cond the differences in density at a single temperature for different compositions. Figure 5 0 shows that as the concentration of lithium is increased the degree of interconnectedness between the lithium sites is visually greater. This is evident of more diffusion occurring at lower temperature s for higher concentration compositions. 131 [Li]=5 [Li]=5.5 [Li]=6 [Li]=6.25 [Li]=6.5 [Li]=6.75 [Li]=7 300 K 600 K 900 K 1200 K Figure 50 : Lithium nuclear density maps for LLZT series along the [100] direction at 300 K, 600 K, 900 K, and 1200 K. Isosurface level 0.3 Å - 3 . 132 One of the major change s among the compositions at 300 K is the shift from more spherical to oblong shaped density regions as lithium - ions shift from tetrahedral to octahedral sites with increasing lithium concentration. The [Li] = 6 composition even shows evidence of long - distance diffusive motion within the volume shown. All compositions except L7LZ show signs of or single jump events with density connecting one tetrahedral and octahedral site. The L7LZ ordering is quite evident at both 300 K and 600 K with identical structures with the only difference being an increase in occupied volume at the higher temperat ure. When the temperature in increased to 600 K for the other compositions, the density maps are much more connected when compared to 300 K. There is now evidence of multiple site jumps for most compositions with the [Li] = 6.5 & 6.75 compositions showing almost complete disorder in agreement with the entropy plots of Figure 49 . The 900 K and 1200 K isosurfaces are completely interconnected for all com positions suggesting that the conductivity weakly composition dependent at high temperatures . Finally, t he density maps at every temper ature confirm the Td - Oh - Td conduction mechanism . This is shown by only having nuclear density within the described sites with diffusion occurring by directly passing through the bottle - neck region connecting the crystallographic sites . The reduced density around the bottle - neck for higher lithium concentration compositions indicates that there are likely more center - pass type events in supporting the cluster analysis and dynamics in Chapter 4. 8.3 Conclusion s T he role of lithium concentration on the structural and dynamic properties in the lithium garnet series was successfully modeled using our co re - only force field. The lattice parameter and 133 phase transformation phenomena were accurately captured as both a function of temperature and lithium showing that the doping of Ta on to Zr sites in L7LZ acts as a means of depressing the phase transition temp erature by introducing lithium disorder and limiting Li Li interactions. The regular solution model was applied to the lithium sublattice with respect to tetrahedral and octahedral sites , showing a trend for decreasing tetrahedral and increasing octahedral occupancy with respect to increasing lithium concentration and temperature. The maximum room temperature conductivity is likely to be observed for the [Li] = 6.5 when the critical lithium concentration is reached . Despite the [Li] = 6.75 composition having a higher self - diffusivity, and conductivity at room temperature , these values are extrapolated from high temperatu re simulations where only cubic LLZT is observed. We believe the tetragonal like structure observed when [Li] = 6.75 indica tes this composition may not be phase stable at room temperature and will likely separate into L6.5LZT and L7LZ phases as suggested by Thompson et al. 149 The excess entropy of the LLZT series was also evaluated, supporting the maximum conductivity at [Li] = 6.5, but will sh ift to [Li] = 6.75 once the transformation to a cubic phase has occurred. Similarly , the excess entropy calculations suggest that the tetragonal to cubic phase transformation for L7LZ is entropy driven, and a first order transformation with respect to the dimensionless heat capacity. 134 C HAPTER 9 : F uture wo rk and c onclusions 9.1 Future work 9.1.1 Exploration in to multivalent doping schemes T he work presented in Chapter 8 suggests that a fundamental way of improving the conductivity in the garnet series is b y maximizing the amount of disorder about the lithium sublatti ce. Having identified that a [Li] = 6.5 is likely the ideal composition for use in devices because it exhibits the highest conductivity, there merits further evaluation into the role of dopants on the lithium conduction pathway by Al 3+ and Ga 3+ substitution , and framework doping onto the La , Ta/Zr and O site s. Based upon the thermodynamic factor and excess entropy results, introducing more lithi um disorder by performing multiple doping strategies with maintaing the critical lithium concentration may result in a more conductive composition. One such theoretical composition would involve doping Al 3+ on to lithium sites, focusing the lithium conduc tion pathway by eliminating oc tahedral and tetrahedral sites. This will likely lower the self - diffusivity, by limiting the possible number of available sites for lithium jumps, but overall may make the collective diffusion higher by forcing more Li Li inte ractions and simultaneous jumps . Subsequence d oping the La site with a 2+ ion and doping Ta/Zr site with a 3+ or 6+ ion, the internal framework of the garnet will become maximally anisotropic. A theoretical composition would look something like Li 6.5 Al 0.16 La 2.35 Y 0.5 Ba 0.15 Zr 1.75 Ta 0.2 Sc 0. 05 O 12 . This composition places multiple atom types on each site effectively changing the bonding character of oxygen within each polyhedral. The introduction of an isotropy within the framework, by initiate continual lithium r earrangements on the lithium the sublattice to obtain a favorable local arrangement. 135 What needs to be further addressed in order to allow such an investigation is the synthesis of such a heavily doped composition. Because so many different elements are in volved, care needs to be taken in order to assure pure phase garnet is obtained through solid - state or sol - gel synthesis. If deemed possible, interatomic potential and DFT model s will need to be created using the techniques used in the previous work . This again introduces problems in the search for an appropriate model , because such a dynamic bonding environment is created, it will be difficult to make it preform well for each pair of atoms. 9.1.2 Electrolyte - electrode interface modeling. W ith the evaluatio n of bulk LLZT reasonably understood, the next major issue to understanding the performance of the garnet as an electrolyte would be to evaluate the electrolyte - electrode interface . To date the formation of a strongly interacting interface remains an elusive problem. It is possible to use spark plasma sintering or sputtering techniques to join the solid electrolyte to an electrode, but the formation of intermediary phases and the qua lity of interface are still relatively unknown. The dynamics of charge transfer across the interface could be evaluated by joining the current model with that of an electrode material to evaluate the conductivity across the region. Impedance spectroscopy h as shown that the density of an electrolyte pellet is critical to having high conductivity with high grain boundary resistances. So, evaluating the quality of interfacial alignment should have similar results and be a major contributor to internal resistan ces. 9.1.3 Evaluation of density correlation functions T he presented work initiates the application of density correlation functions as a means of evaluating the performance of solid electrolytic materials. However, the discrepancy between experimental and simulation results still remain s . In this work we assumed an isotropic response 136 in density fluctuation with respect to our Q vector. It may be that the diffusion in the garnet is not isotropic and the reduction to just using the transverse correlation is insufficient in modeling the system. An exploration of directionally dependent Q vectors therefore merits investigation . A model system could be diffusion in layered materials like lithium cobalt oxide or sodium nickel titanium oxide . Here we would be reduced to two - dimensional planar diffusion, and the role of Q i n directions not parallel to the plane can be investigated. Lastly, because we are no longer assuming isotropic charge density fluctuations, the macroscopic conductivity , , is no longer identical to . 62 The identification of a method that can accurately calculate the conductivity for fa st ion conducting solids is therefore critical to any future work performed within the field of material simulation , and these other systems or methods need be evaluated. 9.2 Conclusions C ompositions within the fast - ionic conducting series of Li 7 - x La 3 Zr 2 - x Ta x O 12 w ere investigated using a combination of molecular dynamics and quasi - elastic neutron scattering techniques. These materials have shown promise as potential solid electrolyte materials in all solid - state lithium - io n batteries , with room temperature ionic conducti vi ties approaching 10 - 3 S/cm . An intimate look at the role of lithium concentration on the observed structural and dynamic properties was evaluated using both classical molecular dynamics and density functio nal theory - based calculations. The analysis techniques employed in this study give rise to best practice s for exploring other fast ion conducting materials with respect to proper model selection and the evaluation of physical properties like conductivity , diffusivity, thermodynamic factor, and entropy. All of which can be easily applied to other materials of interest. 137 First to establish a means of proper model selection, the role the electron exchange - correlation functional (XC) was evaluated on Li 7 La 3 Zr 2 O 12 using density functional theory (DFT) based molecular dynamics simulations. Fourteen different functional forms were evaluated and the prediction of a tetragonal to cubic phase transformation was highly dependent upon what XC functional was chosen. T he functional forms at both small and large density gradient ranges affect ed whether the correct lattice volume or phase behavior was obtained for PBE - like GGA simulations . SOGGA, PBEsol, and PBE2 were the best performing XC functi onal forms generating s tructures comparable to that of experimental of L7LZ . The model garnet materials Li 5 La 3 Ta 2 O 12 (L5LT) and Li 7 La 3 Zr 2 O 12 (L7LZ) were modeled using a core - shell potential model in order to investigate phase change phenomena and the mecha nisms of ionic diffusion with the garnet crystal. The model was able to accurately capture the tetragonal to cubic phase change behavior of L7LZ , predicting a transformation temperature around 900 K , and a consistent cubic phase for L5LT. The lithium - ion dynamic s of these two compositions w ere investigated by observing the role of local clusters on the trajectory of lithium through the bottleneck region connecting the tetrahedral and octahedral sites . Through our observations and analysis by the van Hove correlation function G(r,t), we show that lithium undergoes jump diffusion type dynamics , switching between oscillation within a site and jumping to the nearest neighbor. The modeling also shows that as the concentration of lithium increases more center - pass type diffusion events occur because of the increased interaction of lithium - ions and increased symmetry around tetrahedral sites. From the molecular dynamic simulations on L5LT we were able to calc ulate the intermediate scattering function I(Q,t) and compare the results with experimental dynamic structure factor S(Q, ) measurements using quasi - elastic neutron scattering (QENS) on the 138 BASIS instrument at Oak Ridge National Lab. Our simulations and ex perimental diffusion values are in good agreement with each other predicting a Singwi and S j ö lander (SS) jump type mechanism , quantifying the behavior in G(r,t) performed previously on L7LZ. The relaxed exponential decay of I(Q,t) was fit using the Kohlrau sch - Williams - Watts (KWW) model, with the stretching parameter relating to the degree of interaction between diffusing particles. At low Q, long time and distance, the lithium is weakly interacting ( , while at high Q the values of converge to a plat eau ( < 1) comparable to the distance between two octahedral sites through the tetrahedral. The results of the simulation are in agreement with reports of diffusivity through pulsed gradient spin echo nu clear magnetic resonance ( PGSE - NMR ) and muon sp in re laxation . However, these techniques only resolve residence time values while QENS is able to capture the jump - diffusion distance. Our results suggest more frequent and shorter jumps at high temperature, and slower longer jumps at low temperatures. A core - only classical molecular dynamics model was created to investigate the finite - size effects during simulation and role of lithium concentration on the performance of LLZT as a solid electrolyt e material . The core - only model was able to accurately predict the phase change behavior for L7LZ and lower lithium compositions , suggesting that the introduction of disorder on the lithium sublattice is what introduces the phase change during heating or r eduction in lithium concentration. Based upon the finite - size study, a 3 x 3 x 3, simulation cell is required to obtain convergent behavior in dynamic processes like conductivity, self - diffusivity, Fickian - diffusivity, and thermodynamic factor , so this was determined to be the minimum cell size for evaluation on the effects of lithium concentration on these properties. With respect to lithium concentration, the highest conductivity is obtained when the concentration of lithium is equal to 6.5 p.f.u. , the critical concentration to maintain a single room 139 temperature cubic phase. This maximum shifts to 6.75 as the temperature increases and a single cubic phase is observed for this composition. This is i n agreement with excess entropy calculations showing a identical trend. This suggests that calculating the information entropy of the system may act as a correction factor for non - Arrhenius behavior observed at low temperatures for simulating materials. Th e regular solution model was applied to the lithium sublattice by identifying tetrahedral and octrahedral sites as separate phases. The model predicts a shift from tetrahedral to octahedral site preference as the composition of lithium is increases , minimi zing Li Li interactions, eventually leading to the ordered structure of L7LZ. The thermodynamic factor - 1 , for Li Li interactions tracts closely with that of Li Zr, and opposed to Li Ta pairs, indicating that the framework atoms may play an underlying role in the distribution of lithium leading to a more conductive electrolyte. 140 BIBLIOGRAPHY 141 BIBLIOGRAPHY 1. Wang, Y.; Lai, W., Phase transition in lithium garnet oxide ionic conducto rs Li7La3Zr2O12: the role of Ta substitution and H2O/CO2 exposure. Journal of Power Sources 2015 , 275 , 612 - 620. http://dx.doi.org/10.1016/j.jpowsour.2014.11.062 2. Kamaya, N., et al., A lithium superionic conductor. Nature Materials 2011 , 10 , 682 - 686. http ://dx.doi.org/10.1038/nmat3066 3. Cussen, E. J., The structure of lithium garnets: cation disorder and clustering in a new family of fast Li+ conductors. Chemical Communications 2006 , 412 - 413. http://dx.doi.org/10.1039/b514640b 4. Matsui, M.; Takahashi, K. ; Sakamoto, K.; Hirano, A.; Takeda, Y.; Yamamoto, O.; Imanishi, N., Phase stability of a garnet - type lithium ion conductor Li7La3Zr2O12. Dalton Transactions 2014 , 43 , 1019 - 1024. http://dx.doi.org/10.1039/c3dt52024b 5. Miara, L. J.; Ong, S. P.; Mo, Y. F.; R ichards, W. D.; Park, Y.; Lee, J. M.; Lee, H. S.; Ceder, G., Effect of Rb and Ta Doping on the Ionic Conductivity and Stability of the Garnet Li7+2x - y(La3 - xRbx)(Zr2 - yTay)O - 12 (0 <= x <= 0.375, 0 <= y <= 1) Superionic Conductor: A First Principles Investiga tion. Chemistry of Materials 2013 , 25 , 3048 - 3055. http://dx.doi.org/10.1021/Cm401232r 6. Blomgren, G. E., The Development and Future of Lithium Ion Batteries. Journal of the Electrochemical Society 2017 , 164 , A5019 - A5025. http://dx.doi.org/10.1149/2.025170 1jes 7. Kim, J. G.; Son, B.; Mukherjee, S.; Schuppert, N.; Bates, A.; Kwon, O.; Choi, M. J.; Chung, H. Y.; Park, S., A review of lithium and non - lithium based solid state batteries. Journal of Power Sources 2015 , 282 , 299 - 322. http://dx.doi.org/10.1016/j.j powsour.2015.02.054 8. Aurbach, D.; Talyosef, Y.; Markovsky, B.; Markevich, E.; Zinigrad, E.; Asraf, L.; Gnanaraj, J. S.; Kim, H. J., Design of electrolyte solutions for Li and Li - ion batteries: a review. Electrochim Acta 2004 , 50 , 247 - 254. http://dx.doi.o rg/10.1016/j.electacta.2004.01.090 9. Scrosati, B.; Garche, J., Lithium batteries: Status, prospects and future. Journal of Power Sources 2010 , 195 , 2419 - 2430. http://dx.doi.org/10.1016/j.jpowsour.2009.11.048 10. Scrosati, B.; Hassoun, J.; Sun, Y. K., Lith ium - ion batteries. A look into the future. Energy & Environmental Science 2011 , 4 , 3287 - 3295. http://dx.doi.org/10.1039/c1ee01388b 11. Backstrom, L.; Boldi, P.; Rosa, M.; Ugander, J.; Vigna, S., Four degrees of separation. In Proceedings of the 4th Annual ACM Web Science Conference , ACM: Evanston, Illinois, 2012; pp 33 - 42. 12. Dunn, B.; Kamath, H.; Tarascon, J. M., Electrical Energy Storage for the Grid: A Battery of Choic es. Science 2011 , 334 , 928 - 935. http://dx.doi.org/10.1126/science.1212741 142 13. Goodenough, J. B.; Kim, Y., Challenges for Rechargeable Li Batteries. Chemistry of Materials 2010 , 22 , 587 - 603. http://dx.doi.org/10.1021/cm901452z 14. Bandhauer, T. M.; Garimell a, S.; Fuller, T. F., A Critical Review of Thermal Issues in Lithium - Ion Batteries. Journal of the Electrochemical Society 2011 , 158 , R1 - R25. http://dx.doi.org/10.1149/1.3515880 15. Wang, Q. S.; Ping, P.; Zhao, X. J.; Chu, G. Q.; Sun, J. H.; Chen, C. H., T hermal runaway caused fire and explosion of lithium ion battery. Journal of Power Sources 2012 , 208 , 210 - 224. http://dx.doi.org/10.1016/j.jpowsour.2012.02.038 16. Klenk, M. J.; Lai, W., Finite - size effects on the molecular dynamics simulation of fast - ion c onductors: A case study of lithium garnet oxide Li7La3Zr2O12. Solid State Ionics 2016 , 289 , 143 - 149. http://dx.doi.org/10.1016/j.ssi.2016.03.002 17. Nitta, N.; Wu, F.; Lee, J. T.; Yushin, G., Li - ion battery materials: present and future. Materials Today 20 15 , 18 , 252 - 264. http://dx.doi.org/http://dx.doi.org/10.1016/j.mattod.2014.10.040 18. Xu, K., Nonaqueous liquid electrolytes for lithium - based rechargeable batteries. Chemical Reviews 2004 , 104 , 4303 - 4417. http://dx.doi.org/10.1021/cr030203g 19. Zhang, S. S., A review on electrolyte additives for lithium - ion batteries. Journal of Power Sources 2006 , 162 , 1379 - 1394. http://dx.doi.org/10.1016/j.jpowsour.2006.07.074 20. Edstrom, K.; Gustafsson, T.; Thomas, J. O., The cathode - electrolyte interface in the Li - ion battery. Electrochim Acta 2004 , 50 , 397 - 403. http://dx.doi.org/10.1016/j.electacta.2004.03.049 21. Torabi, F.; Esfahanian, V., Study of Thermal - Runaway in Batteries I. Theoretical Study and Formulation. Journal of the Electrochemical Society 2011 , 158 , A8 50 - A858. http://dx.doi.org/10.1149/1.3592486 22. Oudenhoven, J. F. M.; Baggetto, L.; Notten, P. H. L., All - Solid - State Lithium - Ion Microbatteries: A Review of Various Three - Dimensional Concepts. Advanced Energy Materials 2011 , 1 , 10 - 33. http://dx.doi.org/1 0.1002/aenm.201000002 23. Robinson, A. L.; Janek, J., Solid - state batteries enter EV fray. Mrs Bulletin 2014 , 39 , 1046 - 1047. http://dx.doi.org/10.1557/mrs.2014.285 24. Han, F. D.; Zhu, Y. Z.; He, X. F.; Mo, Y. F.; Wang, C. S., Electrochemical Stability of Li10GeP2S12 and Li7La3Zr2O12 Solid Electrolytes. Advanced Energy Materials 2016 , 6 . http://dx.doi.org/10.1002/aenm.201501590 25. Thangadurai, V.; Kaack, H.; Weppner, W. J. F., Novel fast lithium ion conduction in garnet - type Li5La3M2O12 (M = Nb, Ta). Journ al of the American Ceramic Society 2003 , 86 , 437 - 440. 143 26. Thangadurai, V.; Adams, S.; Weppner, W., Crystal Structure Revision and Identification of Li+ - Ion Migration Pathways in the Garnet - like Li5La3M2O12(M = Nb, Ta) Oxides. Chemistry of Materials 2004 , 1 6 , 2998 - 3006. http://dx.doi.org/10.1021/cm031176d 27. Thangadurai, V.; Weppner, W., Effect of sintering on the ionic conductivity of garnet - related structure Li5La3Nb2O12 and In - and K - doped Li5La3Nb2O12. Journal of Solid State Chemistry 2006 , 179 , 974 - 984 . http://dx.doi.org/10.1016/j.jssc.2005.12.025 28. Murugan, R.; Thangadurai, V.; Weppner, W., Fast lithium ion conduction in garnet - type Li(7)La(3)Zr(2)O(12). Angew Chem Int Ed Engl 2007 , 46 , 7778 - 81. http://dx.doi.org/10.1002/anie.200701144 29. Cussen, E. J., Structure and ionic conductivity in lithium garnets. Journal of Materials Chemistry 2010 , 20 , 5167 - 5173. http://dx.doi.org/10.1039/b925553b 30. Geiger, C. A.; Alekseev, E.; Lazic, B.; Fisch, M.; Armbruster, T.; Langner, R.; Fechtelkord, M.; Kim, N.; P ettke, T.; Weppner, W., Crystal Chemistry and Stability of "Li(7)La(3)Zr(2)O(12)" Garnet: A Fast Lithium - Ion Conductor. Inorganic Chemistry 2011 , 50 , 1089 - 1097. http://dx.doi.org/10.1021/ic101914e 31. Kuhn, A.; Narayanan, S.; Spencer, L.; Goward, G.; Thang adurai, V.; Wilkening, M., Li self - diffusion in garnet - type Li7La3Zr2O12 as probed directly by diffusion - induced Li - 7 spin - lattice relaxation NMR spectroscopy. Physical Review B 2011 , 83 , 094302. http://dx.doi.org/10.1103/PhysRevB.83.094302 32. Teng, S.; T an, J.; Tiwari, A., Recent developments in garnet based solid state electrolytes for thin film batteries. Current Opinion in Solid State & Materials Science 2014 , 18 , 29 - 38. http://dx.doi.org/10.1016/j.cossms.2013.10.002 33. Thangadurai, V.; Narayanan, S.; Pinzaru, D., Garnet - type solid - state fast Li ion conductors for Li batteries: critical review. Chem. Soc. Rev. 2014 , 43 , 4714 - 4727. http://dx.doi.org/10.1039/c4cs00020j 34. Thangadurai, V.; Narayanan, S.; Pinzaru, D., Garnet - type solid - state fast Li ion c onductors for Li batteries: critical review. Chem. Soc. Rev. 2014 , 43 , 4714 - 4727. http://dx.doi.org/10.1039/c4cs00020j 35. Wang, Y. X.; Lai, W., High Ionic Conductivity Lithium Garnet Oxides of Li7 - xLa3Zr2 - xTaxO12 Compositions. Electrochemical and Solid St ate Letters 2012 , 15 , A68 - A71. http://dx.doi.org/10.1149/2.024205esl 36. Xu, M.; Park, M. S.; Lee, J. M.; Kim, T. Y.; Park, Y. S.; Ma, E., Mechanisms of Li+ transport in garnet - type cubic Li3+xLa3M2O12 (M = Te, Nb, Zr). Physical Review B 2012 , 85 , 052301. http://dx.doi.org/10.1103/Physrevb.85.052301 37. Wang, Y. X.; Huq, A.; Lai, W., Insight into lithium distribution in lithium - stuffed garnet oxides through neutron diffraction and atomistic simulation: Li7 - xLa3Zr2 - xTaxO12 (x=0 - 2) series. Solid State Ionics 2014 , 255 , 39 - 49. http://dx.doi.org/10.1016/j.ssi.2013.11.017 144 38. Murugan, R.; Thangadurai, V.; Weppner, W., Fast lithium ion conduction in garnet - type Li7La3Zr2O12. Angew. Chem. - Int. Edit. 2007 , 46 , 7778 - 7781. http://dx.doi.org/10.1002/anie.200701144 39. Awaka, J.; Kijima, N.; Hayakawa, H.; Akimoto, J., Synthesis an d structure analysis of tetragonal Li(7)La(3)Zr(2)O(12) with the garnet - related type structure. Journal of Solid State Chemistry 2009 , 182 , 2046 - 2052. http://dx.doi.org/10.1016/j.jssc.2009.05.020 40. Allen, J. L.; Wolfenstine, J.; Rangasamy, E.; Sakamoto, J., Effect of substitution (Ta, Al, Ga) on the conductivity of Li7La3Zr2O12. Journal of Power Sources 2012 , 206 , 315 - 319. http://dx.doi.org/10.1016/j.jpowsour.2012.01.131 41. Wolfenstine, J.; Ratchford, J.; Rangasamy, E.; Sakamoto, J.; Allen, J. L., Synthe sis and high Li - ion conductivity of Ga - stabilized cubic Li7La3Zr2O12. Materials Chemistry and Physics 2012 , 134 , 571 - 575. http://dx.doi.org/10.1016/j.matchemphys.2012.03.054 42. Thompson, T.; Sharafi, A.; Johannes, M. D.; Huq, A.; Allen, J. L.; Wolfenstine , J.; Sakamoto, J., A Tale of Two Sites: On Defining the Carrier Concentration in Garnet - Based Ionic Conductors for Advanced Li Batteries. Advanced Energy Materials 2015 , 5 . http://dx.doi.org/10.1002/aenm.201500096 43. Larraz, G.; Orera, A.; Sanjuan, M. L. , Cubic phases of garnet - type Li7La3Zr2O12: the role of hydration. Journal of Materials Chemistry A 2013 , 1 , 11419 - 11428. http://dx.doi.org/10.1039/c3ta11996c 44. Bernstein, N.; Johannes, M. D.; Khang, H., Origin of the Structural Phase Transition in Li7La 3Zr2O12. Physical Review Letters 2012 , 109 , 205702. http://dx.doi.org/10.1103/PhysRevLett.109.205702 45. Meier, K.; Laino, T.; Curioni, A., Solid - State Electrolytes: Revealing the Mechanisms of Li - Ion Conduction in Tetragonal and Cubic LLZO by First - Princi ples Calculations. Journal of Physical Chemistry C 2014 , 118 , 6668 - 6679. http://dx.doi.org/10.1021/jp5002463 46. Jalem, R.; Nakayama, M.; Manalastas, W.; Kilner, J. A.; Grimes, R. W.; Kasuga, T.; Kanamura, K., Insights into the Lithium - Ion Conduction Mecha nism of Garnet - Type Cubic Li5La3Ta2O12 by ab - Initio Calculations. Journal of Physical Chemistry C 2015 , 119 , 20783 - 20791. http://dx.doi.org/10.1021/acs.jpcc.5605068 47. Jalem, R.; Yamamoto, Y.; Shiiba, H.; Nakayama, M.; Munakata, H.; Kasuga, T.; Kanamura, K., Concerted Migration Mechanism in the Li Ion Dynamics of Garnet - Type Li7La3Zr2O12. Chemistry of Materials 2013 , 25 , 425 - 430. http://dx.doi.org/10.1021/Cm303542x 48. Jalem, R.; Rushton, M. J. D.; Manalastas, W.; Nakayama, M.; Kasuga, T.; Kilner, J. A.; G rimes, R. W., Effects of Gallium Doping in Garnet - Type Li7La3Zr2O12 Solid Electrolytes. Chemistry of Materials 2015 , 27 , 2821 - 2831. http://dx.doi.org/10.1021/cm5045122 145 49. Rettenwander, D.; Blaha, P.; Laskowski, R.; Schwarz, K.; Bottke, P.; Wilkening, M.; Geiger, C. A.; Amthauer, G., DFT Study of the Role of Al3+ in the Fast Ion - Conductor Li7 - 3xAlx3+La3Zr2O12 Garnet. Chemistry of Materials 2014 , 26 , 2617 - 2623. http://dx.doi.org/10.1021/cm5000999 50. Adams, S.; Rao, R. P., Ion transport and phase transition in Li7 - xLa3(Zr2 - xMx)O - 12 (M = Ta5+, Nb5+, x=0, 0.25). J Mater Chem 2012 , 22 , 1426 - 1434. http://dx.doi.org/10.1039/C1jm14588f 51. Wang, Y. X.; Klenk, M.; Page, K.; Lai, W., Local Structure and Dynamics of Lithium Garnet Ionic Conductors: A Model Material Li 5La3Ta2O12. Chemistry of Materials 2014 , 26 , 5613 - 5624. http://dx.doi.org/10.1021/cm502133c 52. Klenk, M.; Lai, W., Local structure and dynamics of lithium garnet ionic conductors: tetragonal and cubic Li7La3Zr2O12. Physical Chemistry Chemical Physics 2015 , 17 , 8758 - 8768. http://dx.doi.org/10.1039/C4CP05690F 53. Klenk, M. J.; Boeberitz, S. E.; Dai, J.; Jalarvo, N. H.; Peterson, V. K.; Lai, W., Lithium self - diffusion in a model lithium garnet oxide Li5La3Ta2O12: A combined quasi - elastic neutron scattering an d molecular dynamics study. Solid State Ionics 2017 , 312 , 1 - 7. http://dx.doi.org/10.1016/j.ssi.2017.09.022 54. Buschmann, H., et al., Structure and dynamics of the fast lithium ion conductor "Li7La3Zr2O12". Physical Chemistry Chemical Physics 2011 , 13 , 193 78 - 19392. http://dx.doi.org/10.1039/c1cp22108f 55. Duvel, A.; Kuhn, A.; Robben, L.; Wilkening, M.; Heitjans, P., Mechanosynthesis of Solid Electrolytes: Preparation, Characterization, and Li Ion Transport Properties of Garnet - Type Al - Doped Li7La3Zr2O12 Cry stallizing with Cubic Symmetry. Journal of Physical Chemistry C 2012 , 116 , 15192 - 15202. http://dx.doi.org/10.1021/jp301193r 56. Hayamizu, K.; Matsuda, Y.; Matsui, M.; Imanishi, N., Lithium ion diffusion measurements on a garnet - type solid conductor Li6.6La 3Zr1.6Ta0.4O12 by using a pulsed - gradient spin - echo NMR method. Solid State Nuclear Magnetic Resonance 2015 , 70 , 21 - 27. http://dx.doi.org/10.1016/j.ssnmr.2015.05.002 57. Hayamizu, K.; Seki, S.; Haishi, T., Lithium ion micrometer diffusion in a garnet - type cubic Li7La3Zr2O12 (LLZO) studied using Li - 7 NMR spectroscopy. Journal of Chemical Physics 2017 , 146 . http://dx.doi.org/10.1063/1.4973827 58. Wang, D. W.; Zhong, G. M.; Pang, W. K.; Guo, Z. P.; Li, Y. X.; McDonald, M. J.; Fu, R. Q.; Mi, J. X.; Yang, Y., To ward Understanding the Lithium Transport Mechanism in Garnet - type Solid Electrolytes: Li+ Ion Exchanges and Their Mobility at Octahedral/Tetrahedral Sites. Chemistry of Materials 2015 , 27 , 6650 - 6659. http://dx.doi.org/10.1021/acs.chemmater.5b02429 59. Noza ki, H.; Harada, M.; Ohta, S.; Watanabe, I.; Miyake, Y.; Ikedo, Y.; Jalarvo, N. H.; Mamontov, E.; Sugiyama, J., Li diffusive behavior of garnet - type oxides studied by muon - spin 146 relaxation and QENS. Solid State Ionics 2014 , 262 , 585 - 588. http://dx.doi.org/10 .1016/j.ssi.2013.10.014 60. Amores, M.; Ashton, T. E.; Baker, P. J.; Cussen, E. J.; Corr, S. A., Fast microwave - assisted synthesis of Li - stuffed garnets and insights into Li diffusion from muon spin spectroscopy. Journal of Materials Chemistry A 2016 , 4 , 1 729 - 1736. http://dx.doi.org/10.1039/c5ta08107f 61. Han, J. T., et al., Experimental visualization of lithium conduction pathways in garnet - type Li7La3Zr2O12. Chemical Communications 2012 , 48 , 9840 - 9842. http://dx.doi.org/10.1039/c2cc35089k 62. Hansen, J. P .; McDonald, I. R., Theory of Simple Liquids ; Elsevier Ltd, 2013. 63. Mamontov, E.; Herwig, K. W., A time - of - flight backscattering spectrometer at the Spallation Neutron Source, BASIS. Review of Scientific Instruments 2011 , 82 . http://dx.doi.org/10.1063/1. 3626214 64. Arnold, O., et al., Mantid - Data analysis and visualization package for neutron scattering and mu SR experiments. Nuclear Instruments & Methods in Physics Research Section a - Accelerators Spectrometers Detectors and Associated Equipment 2014 , 764 , 156 - 166. http://dx.doi.org/10.1016/j.nima.2014.07.029 65. Azuah, R. T.; Kneller, L. R.; Qiu, Y. M.; Tregenna - Piggott, P. L. W.; Brown, C. M.; C opley, J. R. D.; Dimeo, R. M., DAVE: A Comprehensive Software Suite for the Reduction, Visualization, and Analysis of Low Energy Neutron Spectroscopic Data. Journal of Research of the National Institute of Standards and Technology 2009 , 114 , 341 - 358. http: //dx.doi.org/10.6028/jres.114.025 66. Smith, W.; Forester, T. R., DL_POLY_2.0: A general - purpose parallel molecular dynamics simulation package. Journal of Molecular Graphics 1996 , 14 , 136 - 141. http://dx.doi.org/10.1016/s0263 - 7855(96)00043 - 4 67. Plimpton, S., Fast parallel algorithms for short - range molecular - dynamics. Journal of Computational Physics 1995 , 117 , 1 - 19. http://dx.doi.org/10.1006/jcph.1995.1039 68. Manz, T. A.; Sholl, D. S., Improved Atoms - in - Molecule Charge Partitioning Functional for Simulta neously Reproducing the Electrostatic Potential and Chemical States in Periodic and Nonperiodic Materials. J. Chem. Theory Comput. 2012 , 8 , 2844 - 2867. http://dx.doi.org/10.1021/ct3002199 69. Kresse, G.; Hafner, J., Abinitio molecular - dynamics for liquid - me tals. Physical Review B 1993 , 47 , 558 - 561. http://dx.doi.org/10.1103/PhysRevB.47.558 70. Kresse, G.; Furthmuller, J., Efficiency of ab - initio total energy calculations for metals and semiconductors using a plane - wave basis set. Computational Materials Scie nce 1996 , 6 , 15 - 50. http://dx.doi.org/10.1016/0927 - 0256(96)00008 - 0 147 71. Kresse, G.; Furthmuller, J., Efficient iterative schemes for ab initio total - energy calculations using a plane - wave basis set. Physical Review B 1996 , 54 , 11169 - 11186. http://dx.doi.org /10.1103/PhysRevB.54.11169 72. Kresse, G.; Hafner, J., Ab - initio molecular - dynamics simulation of the liquid - metal amorphous - semiconductor transition in Germanium. Physical Review B 1994 , 49 , 14251 - 14269. http://dx.doi.org/10.1103/PhysRevB.49.14251 73. Kre sse, G.; Joubert, D., From ultrasoft pseudopotentials to the projector augmented - wave method. Physical Review B 1999 , 59 , 1758 - 1775. http://dx.doi.org/10.1103/PhysRevB.59.1758 74. Blochl, P. E., Projector augmented - wave method. Physical Review B 1994 , 50 , 17953 - 17979. http://dx.doi.org/10.1103/PhysRevB.50.17953 75. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized gradient approximation made simple. Physical Review Letters 1996 , 77 , 3865 - 3868. http://dx.doi.org/10.1103/PhysRevLett.77.3865 76. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized gradient approximation made simple (vol 77, pg 3865, 1996). Physical Review Letters 1997 , 78 , 1396 - 1396. http://dx.doi.org/10.1103/PhysRevLett.78.1396 77. Gale, J. D., GULP: A computer program for the symmetry - adap ted simulation of solids. Journal of the Chemical Society - Faraday Transactions 1997 , 93 , 629 - 637. http://dx.doi.org/10.1039/a606455h 78. Gale, J. D.; Rohl, A. L., The General Utility Lattice Program (GULP). Molecular Simulation 2003 , 29 , 291 - 341. http://dx .doi.org/10.1080/0892702031000104887 79. Martyna, G. J.; Tobias, D. J.; Klein, M. L., CONSTANT - PRESSURE MOLECULAR - DYNAMICS ALGORITHMS. Journal of Chemical Physics 1994 , 101 , 4177 - 4189. http://dx.doi.org/10.1063/1.467468 80. Parrinello, M.; Rahman, A., POLY MORPHIC TRANSITIONS IN SINGLE - CRYSTALS - A NEW MOLECULAR - DYNAMICS METHOD. Journal of Applied Physics 1981 , 52 , 7182 - 7190. http://dx.doi.org/10.1063/1.328693 81. Tuckerman, M. E.; Alejandre, J.; Lopez - Rendon, R.; Jochim, A. L.; Martyna, G. J., A Liouville - o perator derived. measure - preserving integrator for molecular dynamics simulations in the isothermal - isobaric ensemble. Journal of Physics a - Mathematical and General 2006 , 39 , 5629 - 5651. http://dx.doi.org/10.1088/0305 - 4470/39/19/s18 82. Shinoda, W.; Shiga, M.; Mikami, M., Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Physical Review B 2004 , 69 . http://dx.doi.org/10.1103/PhysRevB.69.134103 83. Hohenberg, P.; Kohn, W., INHOMOGENEOUS ELECTRON GAS. Physical Review B 1964 , 136 , B864 - +. http://dx.doi.org/10.1103/PhysRev.136.B864 148 84. Kohn, W.; Sham, L. J., SELF - CONSISTENT EQUATIONS INCLUDING EXCHANGE AND CORRELATION EFFECTS. Physical Review 1965 , 140 , 1133 - &. http://dx.doi.org/10.1103/PhysRev.140.A1133 85. Marques, M. A. L.; Oliveira, M. J. T.; Burnus, T., LIBXC: A library of exchange and correlation functionals for density functional theory. Computer Physics Communications 2012 , 183 , 2272 - 2281. http://dx.doi.org/10.1016/j.cpc.2012.05.007 86. Perdew, J. P.; Schmidt, K., Jacob's ladder of density functional approximations for the exchange - correlation energy. In Density Functional Theory and Its Application to Materials , VanDoren, V.; VanAlsenoy, C.; Geerlings, P., Eds. 2001; Vol. 577, pp 1 - 20. 87. Mattsson, A. E.; Armient o, R.; Paier, J.; Kresse, G.; Wills, J. M.; Mattsson, T. R., The AM05 density functional applied to solids. Journal of Chemical Physics 2008 , 128 . http://dx.doi.org/10.1063/1.2835596 88. Constantin, L. A.; Terentjevs, A.; Della Sala, F.; Cortona, P.; Fabia no, E., Semiclassical atom theory applied to solid - state physics. Physical Review B 2016 , 93 . http://dx.doi.org/10.1103/PhysRevB.93.045126 89. Haas, P.; Tran, F.; Blaha, P., Calculation of the lattice constant of solids with semilocal functionals. Physical Review B 2009 , 79 . http://dx.doi.org/10.1103/PhysRevB.79.085104 90. Haas, P.; Tran, F.; Blaha, P., Calculation of the lattice constant of solids with semilocal functionals (vol 79, 085104, 2009). Physical Review B 2009 , 79 . http://dx.doi.org/10.1103/PhysR evB.79.209902 91. Csonka, G. I.; Perdew, J. P.; Ruzsinszky, A.; Philipsen, P. H. T.; Lebegue, S.; Paier, J.; Vydrov, O. A.; Angyan, J. G., Assessing the performance of recent density functionals for bulk solids. Physical Review B 2009 , 79 . http://dx.doi.or g/10.1103/PhysRevB.79.155107 92. Hao, P.; Fang, Y.; Sun, J. W.; Csonka, G. I.; Philipsen, P. H. T.; Perdew, J. P., Lattice constants from semilocal density functionals with zero - point phonon correction. Physical Review B 2012 , 85 . http://dx.doi.org/10.1103/PhysRevB.85.014111 93. Labat , F.; Bremond, E.; Cortona, P.; Adamo, C., Assessing modern GGA functionals for solids. Journal of Molecular Modeling 2013 , 19 , 2791 - 2796. http://dx.doi.org/10.1007/s00894 - 012 - 1646 - 2 94. He, L. H.; Liu, F.; Hautier, G.; Oliveira, M. J. T.; Marques, M. A. L .; Vila, F. D.; Rehr, J. J.; Rignanese, G. M.; Zhou, A. H., Accuracy of generalized gradient approximation functionals for density - functional perturbation theory calculations. Physical Review B 2014 , 89 . http://dx.doi.org/10.1103/PhysRevB.89.064305 95. Ras ander, M.; Moram, M. A., On the accuracy of commonly used density functional approximations in determining the elastic constants of insulators and semiconductors. Journal of Chemical Physics 2015 , 143 . http://dx.doi.org/10.1063/1.4932334 149 96. Tran, F.; Stel zl, J.; Blaha, P., Rungs 1 to 4 of DFT Jacob's ladder: Extensive test on the lattice constant, bulk modulus, and cohesive energy of solids. Journal of Chemical Physics 2016 , 144 . http://dx.doi.org/10.1063/1.4948636 97. Wang, Y. X.; Lai, W., Phase transitio n in lithium garnet oxide ionic conductors Li7La3Zr2O12: The role of Ta substitution and H2O/CO2 exposure. Journal of Power Sources 2015 , 275 , 612 - 620. http://dx.doi.org/10.1016/j.jpowsour.2014.11.062 98. Matsuda, Y.; Sakamoto, K.; Matsui, M.; Yamamoto, O. ; Takeda, Y.; Imanishi, N., Phase formation of a garnet - type lithium - ion conductor Li - 7 ( - ) 3xAlxLa3Zr2O12. Solid State Ionics 2015 , 277 , 23 - 29. http://dx.doi.org/10.1016/j.ssi.2015.04.011 99. Kang, S. G.; Sholl, D. S., First - Principles Study of Chemical S tability of the Lithium Oxide Garnets Li7La3M2O12 (M = Zr, Sn, or Hf). Journal of Physical Chemistry C 2014 , 118 , 17402 - 17406. http://dx.doi.org/10.1021/jp504314w 100. Haas, P.; Tran, F.; Blaha, P.; Pedroza, L. S.; da Silva, A. J. R.; Odashima, M. M.; Cape lle, K., Systematic investigation of a family of gradient - dependent functionals for solids. Physical Review B 2010 , 81 . http://dx.doi.org/10.1103/PhysRevB.81.125136 101. Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constant in, L. A.; Zhou, X. L.; Burke, K., Restoring the density - gradient expansion for exchange in solids and surfaces. Physical Review Letters 2008 , 100 . http://dx.doi.org/10.1103/PhysRevLett.100.136406 102. Hammer, B.; Hansen, L. B.; Norskov, J. K., Improved ad sorption energetics within density - functional theory using revised Perdew - Burke - Ernzerhof functionals. Physical Review B 1999 , 59 , 7413 - 7421. http://dx.doi.org/10.1103/PhysRevB.59.7413 103. Zhao, Y.; Truhlar, D. G., Construction of a generalized gradient a pproximation by restoring the density - gradient expansion and enforcing a tight Lieb - Oxford bound. Journal of Chemical Physics 2008 , 128 . http://dx.doi.org/10.1063/1.2912068 104. Ruzsinszky, A.; Csonka, G. I.; Scuseria, G. E., Regularized Gradient Expansion for Atoms, Molecules, and Solids. J. Chem. Theory Comput. 2009 , 5 , 763 - 769. http://dx.doi.org/10.1021/ct8005369 105. Wu, Z. G.; Cohen, R. E., More accurate generalized gradient approximation for solids. Physical Review B 2006 , 73 . http://dx.doi.org/10.110 3/PhysRevB.73.235116 106. Sarmiento - Perez, R.; Botti, S.; Marques, M. A. L., Optimized Exchange and Correlation Semilocal Functional for the Calculation of Energies of Formation. J. Chem. Theory Comput. 2015 , 11 , 3844 - 3850. http://dx.doi.org/10.1021/acs.jc tc.5b00529 107. Vela, A.; Pacheco - Kato, J. C.; Gazquez, J. L.; del Campo, J. M.; Trickey, S. B., Improved constraint satisfaction in a simple generalized gradient approximation exchange functional. Journal of Chemical Physics 2012 , 136 . http://dx.doi.org/1 0.1063/1.3701132 150 108. VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J., QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Computer Physics Communications 2005 , 1 67 , 103 - 128. http://dx.doi.org/10.1016/j.cpc.2004.12.014 109. Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J., CP2K: atomistic simulations of condensed matter systems. Wiley Interdisciplinary Reviews - Computational Molecular Science 2014 , 4 , 15 - 2 5. http://dx.doi.org/10.1002/wcms.1159 110. Goedecker, S.; Teter, M.; Hutter, J., Separable dual - space Gaussian pseudopotentials. Physical Review B 1996 , 54 , 1703 - 1710. http://dx.doi.org/10.1103/PhysRevB.54.1703 111. Krack, M., Pseudopotentials for H to Kr optimized for gradient - corrected exchange - correlation functionals. Theoretical Chemistry Accounts 2005 , 114 , 145 - 152. http://dx.doi.org/10.1007/s00214 - 005 - 0655 - y 112. VandeVondele, J.; Hutter, J., Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. Journal of Chemical Physics 2007 , 127 . http://dx.doi.org/10.1063/1.2770708 113. Bussi, G.; Donadio, D.; Parrinello, M., Canonical sampling through velocity rescaling. Journal of Chemical Physics 2007 , 126 . http://dx.do i.org/10.1063/1.2408420 114. Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L., Explicit reversible integrators for extended systems dynamics. Molecular Physics 1996 , 87 , 1117 - 1157. http://dx.doi.org/10.1080/00268979650027054 115. Perdew, J. P. ; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C., ATOMS, MOLECULES, SOLIDS, AND SURFACES - APPLICATIONS OF THE GENERALIZED GRADIENT APPROXIMATION FOR EXCHANGE AND CORRELATION. Physical Review B 1992 , 46 , 6671 - 6687 . http://dx.doi.org/10.1103/PhysRevB.46.6671 116. Odashima, M. M.; Capelle, K., How tight is the Lieb - Oxford bound? Journal of Chemical Physics 2007 , 127 . http://dx.doi.org/10.1063/1.2759202 117. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H., A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT - D) for the 94 elements H - Pu. Journal of Chemical Physics 2010 , 132 . http://dx.doi.org/10.1063/1.3382344 118. Grimme, S.; Ehrlich, S.; Goerigk, L., Effect of the Dampin g Function in Dispersion Corrected Density Functional Theory. Journal of Computational Chemistry 2011 , 32 , 1456 - 1465. http://dx.doi.org/10.1002/jcc.21759 119. Jobic, H.; Theodorou, D. N., Quasi - elastic neutron scattering and molecular dynamics simulation as complementary techniques for studying diffusion in zeolites. Microporous and Mesoporous Materials 2007 , 102 , 21 - 50. http://dx.doi.org/10.1016/j.micromeso.2006.12.034 151 120. Kubo, R ., STATISTICAL - MECHANICAL THEORY OF IRREVERSIBLE PROCESSES .1. GENERAL THEORY AND SIMPLE APPLICATIONS TO MAGNETIC AND CONDUCTION PROBLEMS. Journal of the Physical Society of Japan 1957 , 12 , 570 - 586. http://dx.doi.org/10.1143/jpsj.12.570 121. Chudley, C. T. ; Elliott, R. J., NEUTRON SCATTERING FROM A LIQUID ON A JUMP DIFFUSION MODEL. Proceedings of the Physical Society of London 1961 , 77 , 353 - &. http://dx.doi.org/10.1088/0370 - 1328/77/2/319 122. Hall, P. L.; Ross, D. K., INCOHERENT NEUTRON - SCATTERING FUNCTIONS FOR RANDOM JUMP DIFFUSION IN BOUNDED AND INFINITE MEDIA. Molecular Physics 1981 , 42 , 673 - 682. http://dx.doi.org/10.1080/00268978100100521 123. Singwi, K. S.; Sjolander, A., RESONANCE ABSORPTION OF NUCLEAR GAMMA - RAYS AND THE DYNAMICS OF ATOMIC MOTIONS. Phy sical Review 1960 , 120 , 1093 - 1102. http://dx.doi.org/10.1103/PhysRev.120.1093 124. Fuchs, M., THE KOHLRAUSCH LAW AS A LIMIT SOLUTION TO MODE - COUPLING EQUATIONS. Journal of Non - Crystalline Solids 1994 , 172 , 241 - 247. http://dx.doi.org/10.1016/0022 - 3093(94)90 442 - 1 125. Siepmann, J. I.; McDonald, I. R.; Frenkel, D., FINITE - SIZE CORRECTIONS TO THE CHEMICAL - POTENTIAL. Journal of Physics - Condensed Matter 1992 , 4 , 679 - 691. http://dx.doi.org/10.1088/0953 - 8984/4/3/009 126. Schnell, S. K.; Liu, X.; Simon, J. M.; Bardo w, A.; Bedeaux, D.; Vlugt, T. J. H.; Kjelstrup, S., Calculating Thermodynamic Properties from Fluctuations at Small Scales. Journal of Physical Chemistry B 2011 , 115 , 10911 - 10918. http://dx.doi.org/10.1021/jp204347p 127. Schnell, S. K.; Vlugt, T. J. H.; Si mon, J. - M.; Bedeaux, D.; Kjelstrup, S., Thermodynamics of small systems embedded in a reservoir: a detailed analysis of finite size effects. Molecular Physics 2012 , 110 , 1069 - 1079. http://dx.doi.org/10.1080/00268976.2011.637524 128. Compaan, K.; Haven, Y., CORRELATION FACTORS FOR DIFFUSION IN SOLIDS. Trans. Faraday Soc. 1956 , 52 , 786 - 801. http://dx.doi.org/10.1039/tf9565200786 129. Kirkwood, J. G.; Buff, F. P., THE STATISTICAL MECHANICAL THEORY OF SOLUTIONS .1. Journal of Chemical Physics 1951 , 19 , 774 - 777. http://dx.doi.org/10.1063/1.1748352 130. Kruger, P.; Schnell, S. K.; Bedeaux, D.; Kjelstrup, S.; Vlugt, T. J. H.; Simon, J. M., Kirkwood - Buff Integrals for Finite Volumes. Journal of Physical Chemistry Letters 2013 , 4 , 235 - 238. http://dx.doi.org/10.1021/j z301992u 131. Giorgini, S.; Pitaevskii, L. P.; Stringari, S., Anomalous fluctuations of the condensate in interacting Bose gases. Physical Review Letters 1998 , 80 , 5040 - 5043. http://dx.doi.org/10.1103/PhysRevLett.80.5040 152 132. Shannon, C. E., A MATHEMATICAL THEORY OF COMMUNICATION. Bell System Technical Journal 1948 , 27 , 623 - 656. http://dx.doi.org/10.1002/j.1538 - 7305.1948.tb00917.x 133. Yu, H. S.; He, X.; Truhlar, D. G., MN15 - L: A New Local Exchange - Correlation Functional for Kohn - Sham Density Functional The ory with Broad Accuracy for Atoms, Molecules, and Solids. J. Chem. Theory Comput. 2016 , 12 , 1280 - 1293. http://dx.doi.org/10.1021/acs.jctc.5b01082 134. Sun, J. W., et al., Accurate first - principles structures and energies of diversely bonded systems from an efficient density functional. Nature Chemistry 2016 , 8 , 831 - 836. http://dx.doi.org/10.1038/nchem.2535 135. Keen, D. A., Disordering phenomena in superionic conductors. Journal of Physics - Condensed Matter 2002 , 14 , R819 - R857. http://dx.doi.org/10.1088/0953 - 8984/14/32/201 136. Katharina Meier, T. L., and Alessandro Curioni, Solid - State Electrolytes: Revealing the Mechanisms of Li - Ion Conduction in Tetragonal and Cubic LLZO by First - Principles Calculations. J. Phys. Chem. C 2014 , 118 , 6668 - 6679. http://dx.doi .org/10.1021/jp5002463 137. Xie, H.; Alonso, J. A.; Li, Y. T.; Fernandez - Diaz, M. T.; Goodenough, J. B., Lithium Distribution in Aluminum - Free Cubic Li7La3Zr2O12. Chemistry of Materials 2011 , 23 , 3587 - 3589. http://dx.doi.org/10.1021/cm201671k 138. Singwi, K. S.; Sjolander, A., DIFFUSIVE MOTIONS IN WATER AND COLD NEUTRON SCATTERING. Physical Review 1960 , 119 , 863 - 871. http://dx.doi.org/10.1103/PhysRev.119.863 139. Williams, G.; Watts, D. C., NON - SYMMETRICAL DIELECTRIC RELAXATION BEHAVIOUR ARISING FROM A SIMP LE EMPIRICAL DECAY FUNCTION. Trans. Faraday Soc. 1970 , 66 , 80 - +. http://dx.doi.org/10.1039/tf9706600080 140. Trachenko, K.; Brazhkin, V. V., Collective modes and thermodynamics of the liquid state. Reports on Progress in Physics 2016 , 79 . http://dx.doi.org /10.1088/0034 - 4885/79/1/016502 141. Ngai, K. L., Short - time and long - time relaxation dynamics of glass - forming substances: a coupling model perspective. Journal of Physics - Condensed Matter 2000 , 12 , 6437 - 6451. http://dx.doi.org/10.1088/0953 - 8984/12/29/316 142. Sciortino, F.; Gallo, P.; Tartaglia, P.; Chen, S. H., Supercooled water and the kinetic glass transition. Physical Review E 1996 , 54 , 6331 - 6343. http://dx.doi.org/10.1103/PhysRevE.54.6331 143. Yeh, I. C.; Hummer, G., System - size dependence of diffusio n coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions. Journal of Physical Chemistry B 2004 , 108 , 15873 - 15879. http://dx.doi.org/10.1021/jp0477147 144. Zhang, Y. G.; Guo, G. J.; Refson, K.; Zhao, Y. J., Finite - size effect at both high and low temperatures in molecular dynamics calculations of the self - diffusion coefficient and viscosity of 153 liquid silica. Journal of Physics - Condensed Matter 2004 , 16 , 9127 - 9135. http://dx.doi.org/10.1088/0953 - 8984/16/50/003 145. Heyes, D. M.; Cass, M. J.; Powles, J. G.; Evans, W. A. B., Self - diffusion coefficient of the hard - sphere fluid: System size dependence and empirical correlations. Journal of Physical Chemistry B 2007 , 111 , 1455 - 1464. http://dx.doi.org/10.1021/jp067373s 146. Dunweg, B.; Kremer, K., MOLECULAR - DYNAMICS SIMULATION OF A POLYMER - CHAIN IN SOLUTION. Journal of Chemical Physics 1993 , 99 , 6983 - 6997. 147. Villamaina, D.; Trizac, E., Thinking outside the box: fluctuations and finite size effects. European Journal of Physics 2014 , 35 . http://dx.doi.org/10.1088/0143 - 0807/35/3/035011 148. Klenk, M.; Lai, W., Local structure and dynamics of lithium garnet ionic conductors: tetragonal and cubic Li7La3Zr2O7. Phys Chem Chem Phys 2015 , 17 , 8758 - 68. http://dx.doi.org/10.1039/ c4cp05690f 149. Thompson, T.; Wolfenstine, J.; Allen, J. L.; Johannes, M.; Huq, A.; David, I. N.; Sakamoto, J., Tetragonal vs. cubic phase stability in Al - free Ta doped Li7La3Zr2O12 (LLZO). Journal of Materials Chemistry A 2014 , 2 , 13431 - 13436. http://dx .doi.org/10.1039/c4ta02099e 150. Dyre, J. C.; Hechsher, T.; Niss, K., A brief critique of the Adam - Gibbs entropy model. Journal of Non - Crystalline Solids 2009 , 355 , 624 - 627. http://dx.doi.org/10.1016/j.jnoncrysol.2009.01.039