STRUCTURE AND PROPERTIES STUDY ON ENERGY MATERIALS: THERMOELECTRIC MATERIAL TETRAHEDRITE AND LITHIUM ION CONDUCTOR LiPON By Junchao Li A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of M aterials Science and Engineering Doctor of Philosophy 20 20 A BSTRACT STRUCTURE AND PROPERTIES STUDY ON ENERGY MATERIALS: THERMOELECTRIC MATERIAL TETRAHEDRITE AND LITHIUM ION CONDUCTOR LiPON B y Junchao Li Development of efficient energy materials is critic al in order to ease the energy demand and reduce our dependence of fos sil fuel . Thermoelectric materials are promising due to their capability of generating electrical power by recovering waste heat. The performance of thermoelectric materials is quantifie d by a dimensionless figure of merit zT , which depends on their properties such as electrical conductivity, Seebeck coefficient and thermal conductivity. Tetrahedrites, a copper antimony sulfosalt mineral, typified by Cu 12 - x M x Sb 4 S 13 , where M is a transitio n metal element such as Ni, Zn, Fe or Mn, have great potential f or thermoelectric application due to their relatively high zT (close to 1 at 700 K), earth - abundance, environmental friendliness, favorable electrical properties, and most importantly intrinsi c low lattice thermal conductivity (less than 1 W m - 1 K - 1 ) in wide temperature . In addition to energy recovery, reliable energy storage devices are also emerging to relieve the energy demand and improve the efficiency of consuming energy resources. Lithiu m - ion batteries are known to be reliable and successful electr ochemical energy storage devices and appliable in various aspects , including laptops, smartphones and electrical vehicles. Lithium phosphorous oxynitride ( LiPON ) are widely used as thin - film sol id - state electrolytes in Li - ion battery , which is the only dem onstrated solid - state electrolyte that is quite stable in direct contact with Li metal at potentials from 0 - 5 V . However, the structure of LiPON , the effects of N doping, and the origin of its g ood electrochemical stability remains inconclusive. In this thesis, reliable modeling techniques accompanied with experimental tools, are applied to study the thermoelectric material tetrahedrite and the ionic conductor LiPON , in order to study their struc tural and dynamical properties. Accurate and efficient d ensity - functional theory (DFT) and density - functional tight - binding (DFTB) met hods, combined with molecular dynamics (MD) simulation s are utilized in order to investigate the structures and properties of these energy materials. The incoherent and coherent atomic dynamics study of tetrahedrite Cu 10.5 NiZn 0.5 Sb 4 S 13 provides the origin of softening upon cooling by investigate the motion of Cu 12e at different temperatures. The dynamic structure factors in t he longitudinal and transverse direction will also be discussed. The Cu movement of Cu - rich tetrahedrite Cu 1 4 Sb 4 S 13 is revealed by Cu self - diffusivity , nuclear density map and . Moreover, we investig ate the effect of simulation c ell size and basis sets on the DFT - based MD simulation results using tetrahedrite Cu 10 Zn 2 Sb 4 S 13 thermoelectric as a model material , showing the advantage of larger cell by accessing smaller Q range. I n addition, the low - temperature structural properties of Cu 1 2 Sb 4 S 13 is measured by neutron diffraction, which indicates that no cubic to tetragonal transition occurs at metal - semiconductor transition (MST) temperature . Thermoelectric properties such as See beck coefficient, electrical resistivity and electrical thermal conductivity will also be investigated. DFTB method is implemented to study the structure and transport properties of Li 3 PO 4 and LiPON , while the exploration of N doping effect is included. Lastly, the LiPON /Li interphase will be revealed in order to study the origin of electrochemical stability. iv Dedicated to Families and friends v ACKNOWLEDGEMENTS First and foremost, I would like to thank my advisor Dr. Wei Lai for his generous support throughout my path on pursuing doctorate degree. Over these ye ars, He is not only my research advisor but also a spiritual mentor, with his valuable characteristics of hardworking, patience and intelligence. The most important thing s I lea r ned from him is self - motivation and the capability of critical thinking. Th e weekly research update and presentations in our group help me a lot in t erms of oral presentation skill, as well as confidence. I want to show my great appreciation to Dr. Lai again for everything he did these years. I feel so lucky joining our group. I also want to thank my committee members, Dr. Donald Morelli , Dr. Yue Qi an d Dr. helped me a lot and inspired my research path. Dr. Ke, thank you for your valua ble advice during our collaboration on thermoelectric materials and my com prehensive exam. Thank you all again for being my committee members. I would like to ack now ledge my former and current colleagues, Rengarajan Shanmugam, Matt h ew Klenk, Qian Chen, Jin Dai and Yue Jiang. Although it was only half a year working with Raj, he showed a good model as a hardworking PhD candidate. Matt helped me a lot with the computational setup and codes writing when I first joined the group. The discussion between us alway s kept me motivated and brought me new ideas. Qian and Jin, they are patie nt, unselfish and optimistic. Although we are working on different material systems, the works you did inspired me on my own research. I wish you will have successful career at MSU a nd in the future. Yue made a great effort on his experimental work and he is always accommodating. I hope you enjoy your time in our group. vi thank you for everythin g you did throughout the 28 years. Every time I felt frustrated, you alway s My grandmother, thank you for raising me up and kept me accompanied in my childhood. My aunts and uncles, thank you for your concerns and making me feel warm e very time I come back home. You are all the greatest families. I would also like to acknowledge my friends at MSU. We spend many wonderful times together , playing basketball, swimmin g, traveling, etc. The great times we had, the difficulties we got over, and worthy suggestion s you gave, th ose all provide d me motivation on being a better person. vii TABLE OF CONTENTS LIST OF T A BLES ................................ ................................ ................................ .......................... x LIST OF FIGURES ................................ ................................ ................................ ....................... xi CHAPTER 1 : Introd uction and motivation ................................ ................................ ................... 1 1 .1 E nergy crisis ................................ ................................ ................................ .......................... 1 1.2 Thermoelectric materials and tetrahedrite Cu 12 Sb 4 S 13 ................................ .......................... 3 1.2.1 Thermoelectric effect and performance of thermoe lectri c materials .............................. 3 1.2.2 Applications of thermoelectric materials ................................ ................................ ........ 6 1.2.3 Tetrahedrite Cu 12 Sb 4 S 13 ................................ ................................ ................................ .. 7 1.3 Lithium - ion battery and electrolyte material LiPON ................................ .......................... 11 1.3.1 Introduction of lithium - ion battery ................................ ................................ ............... 11 1.3.2 Ele ctrolyte material LiPON ................................ ................................ .......................... 15 1.4 Modeling techniques on structural and dynamical properties of energy materials ............. 18 1.4.1 In troduction to DFT ................................ ................................ ................................ ...... 18 1.4.2 Introduction to DFTB ................................ ................................ ................................ ... 22 1.4.3 Other computational techniques ................................ ................................ ................... 26 1.5 Dissertation outline ................................ ................................ ................................ ............. 27 CHAPTER 2 : Pristine tetrahedrite Cu 12 Sb 4 S 13 : neutron diffraction and computational study on low - temperature properties ................................ ................................ ................................ ........... 31 2.1 Introduction ................................ ................................ ................................ ......................... 31 2.2 Computational and experimental details ................................ ................................ ............. 32 2.3 Results and Discussion ................................ ................................ ................................ ........ 33 2.3.1 Low temperature structure of Cu 12 Sb 4 S 13 ................................ ................................ ..... 33 2.3.2 Sb - Cu 12e - Sb interaction ................................ ................................ .............................. 38 2.3.3 Vibrational Density of States ................................ ................................ ........................ 40 2.3.4 Thermoelectric properties ................................ ................................ ............................. 41 2.4 Conclusion ................................ ................................ ................................ ........................... 45 CHAPTER 3 : Incoherent and coherent atomic dynamics of Cu 10.5 NiZn 0.5 Sb 4 S 13 ...................... 47 3.1 Introduction ................................ ................................ ................................ ......................... 47 3.2 Computational and e xperimental details ................................ ................................ ............. 47 3.3 Results and discussion ................................ ................................ ................................ ......... 48 3.3.1 Incoherent atomic dynamics of Cu 10.5 NiZn 0.5 Sb 4 S 13 ................................ .................... 48 3.3.2 Coherent atomic dynamics of Cu 10.5 NiZn 0.5 Sb 4 S 13 ................................ ...................... 52 3.4 Summary ................................ ................................ ................................ ............................. 53 CHAPTER 4 : Mobile Cu m ovement in Cu - rich tetrahedrite Cu 14 Sb 4 S 13 ................................ .... 55 4.1 Introduction ................................ ................................ ................................ ......................... 55 4.2 Computationa l details ................................ ................................ ................................ .......... 56 viii 4.3 Results and discussion ................................ ................................ ................................ ......... 57 4.3.1 Lattice parameters ................................ ................................ ................................ ......... 57 4.3.2 Incoherent density correlation and Cu self - diffu sivity ................................ ................. 58 4.3.3 Cu nuclear density map ................................ ................................ ................................ 60 4.3.4 Energy barriers ................................ ................................ ................................ ............. 62 4.4 Su mmary ................................ ................................ ................................ ............................. 63 CHAPTER 5 : Density - functional theory based molecular dynamics simulation of tetrahedrite thermoelectrics: effect of cell size and basis sets ................................ ................................ .......... 65 5.1 Introduction ................................ ................................ ................................ ......................... 65 5.2 Computational and experimental details ................................ ................................ ............. 66 5.3 Results and discussion ................................ ................................ ................................ ......... 67 5.3.1 Average structure ................................ ................................ ................................ .......... 67 5.3.2 Latt ice parameters ................................ ................................ ................................ ......... 68 5.3.3 Vibrational density of states ................................ ................................ ......................... 69 5.3.4 Phonon dispersion ................................ ................................ ................................ ......... 71 5.3.5 EXAF S spectra ................................ ................................ ................................ ............. 72 5.3.6 2x2x2 AO Simul ation ................................ ................................ ................................ ... 73 5.4 Summary ................................ ................................ ................................ ............................. 76 CHAPTER 6 : Structure and ionic conduction study on Li 3 PO 4 and LiPON with density functional tight binding (DF TB) method ................................ ................................ ...................... 78 6.1 Introduction ................................ ................................ ................................ ......................... 78 6.2 Computational detail s ................................ ................................ ................................ .......... 79 6.2.1 Gene ration of c - Li 3 PO 4 and c - LiPON models ................................ .............................. 79 6.2.2 Generation of a - Li 3 PO 4 and a - LiPON models ................................ .............................. 79 6.3 Results and discussion ................................ ................................ ................................ ......... 81 6.3.1 Lattice parameters ................................ ................................ ................................ ......... 81 6.3.2 Exploration of LiPON structure ................................ ................................ ................... 81 6.3.3 Self - diffus ivity of different atomic groups ................................ ................................ ... 85 6.3.4 Ionic conductivity ................................ ................................ ................................ ......... 90 6.4 Sum mary ................................ ................................ ................................ ............................. 92 CHAPTER 7 : Preliminary study on LiPON/Li interphase ................................ .......................... 94 7.1 Introduction ................................ ................................ ................................ ......................... 94 7.2 Computational details ................................ ................................ ................................ .......... 94 7.3 Results and discussion ................................ ................................ ................................ ......... 95 7.3.1 Average structure ................................ ................................ ................................ .......... 95 7.3.2 Reac tion at the interface ................................ ................................ ............................... 96 7.3.3 Li charge distribution ................................ ................................ ................................ .... 97 7.3.4 Projected density of states ................................ ................................ ............................ 99 CHAPTER 8 : Future work and conclusions ................................ ................................ .............. 101 8.1 Future work ................................ ................................ ................................ ....................... 101 8.1. 1 Further investigation on thermoelectric material tetrahedrite ................................ .... 101 8.1.2 Understanding the origin of LiPON/Li interphase ................................ ..................... 101 ix 8.1.3 Exploration in cathode materials ................................ ................................ ................ 102 8.2 Conclusions ................................ ................................ ................................ ....................... 102 BIBLIOGRAPHY ................................ ................................ ................................ ....................... 107 x L IST OF TABLES Table 2.1 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 20 K. ............................ 37 Table 2.2 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 50 K. ............................ 37 Tabl e 2.3 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 87.5 K. ......................... 37 Table 2.4 : Atomic parameters in t he crystal structure of Cu 12 Sb 4 S 13 at 300 K. .......................... 37 Table 2.5 : Number of rattling and locking Cu 12e, Cu 12e - Sb bond order and bond distance for Cu 12 Sb 4 S 13 at different temperatures. ................................ ................................ ........................... 39 Table 6.1: Calculated lattice pa rameters (a, b, c) of c - Li 3 PO 4 and c - LiPON at 300 K comparing to experiments. ................................ ................................ ................................ .............................. 81 xi LIST OF FIGURES Figure 1.1: Global energy consumption from 1993 to 2018 wit h the unit of mtoe (million tonnes oil equivalent) 1 . ................................ ................................ ................................ ............................... 1 Figure 1.2: Global energy - related carbon dioxide emiss ions by source, 1990 - 2018 2 . ................... 2 Figure 1.3 : Schematic of thermoelectric device showing the direction of charge flow on both cooling and power generation 11 . ................................ ................................ ................................ ..... 4 Figure 1.4: (a) Optimizing through carrier concentration tu ning. (b) Benefits of reducing . 11 ................................ ................................ ................................ ................................ ......................... 5 Figure 1.5: Crystal structure of tetrahedrite Cu 12 Sb 4 S 13 with two distinct Cu sites. ..................... 8 Figur e 1.6: Electronic band structure and density of states of (a) Cu 12 Sb 4 S 13 and (b) Cu 10 Zn 2 Sb 4 S 13 35 . ................................ ................................ ................................ ........................... 10 Figure 1.7: Schematic of working principles of Li x C 6 /L i 1 - x CoO 2 Li - ion battery 10 . ...................... 12 Figure 1.8: Current electrode materials road map 67 . ................................ ................................ ..... 13 Figure 1.9: Temperature dependence of ionic con ductivity of various electrolyte materials 72 . ... 14 Figure 1.10: Composition diagram with starting materials of LiO 1/2 , LiN 1/3 , PO 5/2 and PN 5/3 . Natural and synthetic crystalline materials are labe led as dark (bl ue) circles, thin - film compo sitions as light (turquoise) circles, and examples of stable and meta - stables nitride phosphate materials as black square 83 . ................................ ................................ .......................... 16 Figure 1.11: Current - po tential curve of a LiPON/Pt half - cell with r espect to Li/Li + reference 91 . 17 Fi gure 1.12: Popular computational techniques for different length and time scales. .................. 18 Figure 1.13: T he number of DFT citations from 19 95 to 2013 102 . ................................ ............... 21 Figure 2.1 : Neutron diffraction data for Cu 12 Sb 4 S 13 at 20, 50, 87.5 and 300 K, synchrontron x - r ay data by Nasonova el al. are shown for compar ison. No peak splitting can be identified below MST temperature. ................................ ................................ ................................ ......................... 34 Figure 2.2 : Neutron diffraction data collected for tetrahedrite sample at 87.5 K using (a, b) Cu 12 Sb 4 S 13 and (c,d) Cu 13.2 Sb 4 S 1 3 as model structures. Block 1 contains 80% of the data and block 2 has the rest 20%. ................................ ................................ ................................ .............. 35 Figure 2.3 : Neutron diffraction data collec ted for tetra hedrite sample at (a, b) 20 K, (c, d) at 50 K and (e, f) at 300 K using Cu 12 Sb 4 S 13 structure as model. Block 1 contains 80% of the data and block 2 has the rest 20%. ................................ ................................ ................................ .............. 36 xii Figure 2.4 : Lattice para meters of Cu 12 Sb 4 S 13 obtained by DFT - MD NPT simula tion and neutron diffraction refinement. Lattice parameters proposed by Nasonova using synchrontron x - ray are shown here for comparison. ................................ ................................ ................................ .......... 38 Figu re 2.5 : Schematic of Cu - Sb interaction at 20 K and 300 K (a) and potential energy landscape based on Sb - Cu12e - Sb cluster at (b) 20 K and (c) 300 K. ................................ ............................ 40 Figure 2.6 : Partial vibration density of s tates (VDOS) of each atomic group for Cu 12 Sb 4 S 13 at different temperatures. ................................ ................................ ................................ .................. 41 Figure 2.7 : DOS calculation for Cu 12 Sb 4 S 13 and partial contribution from each atomic group at different temperatur es. Fermi energy is marked by dashed line. ................................ .................. 42 Figure 2.8 : Temperature dependencies of calculated (a) electrical resistivity, (b) Seebeck coefficient and (c) electrical thermal c onductivity for Cu 12 Sb 4 S 13 . Experimental data from Nasonova el al. 136 are shown here for comparison. ................................ ................................ ...... 44 Figure 3.1 : (a)Comparison of VDOS from MD simulation and INS experiment. (b) Partial VDOS for each atom gr oup and total VDOS from MD simulation at 300K. ............................... 49 Figure 3.2 : Partial VDOS of Cu 12e atoms at different temperatures. Shape was fitted to sum of two Lorentzian functions. ................................ ................................ ................................ ............. 50 Figure 3.3 : Pote ntial energy landscape based on Sb - Cu 12e - Sb cluster at 50 K (a) and 300 K (b), and correlated VDOS spectra of Cu 12e at 50 K (c) and 300 K (d). Cu 12e - Sb interaction are shown in (e). ................................ ................................ ................................ ................................ .. 51 Figure 3.4 : Calculated dynamic structure factors based on the momentum correlation function (a) longitudinal and (b) transverse direction. ................................ ................................ ..................... 52 Figu re 4.1: Initial crystal structure of Cu - rich tetrahedrite C u 14 Sb 4 S 13 with three distinct Cu sites. Cu24g atoms are displayed as black. ................................ ................................ ............................ 56 Figure 4.2: Cu24g local environm ent for better visualiza tion. ................................ ..................... 56 Figure 4.3 : Lattice parameters (a, b, c) as function of time in NPT simulation. .......................... 57 Figure 4.4 : Incoherent density correlati on for (a) diffusing Cu atoms and (b) vibrating Cu atoms for t hree small Q values. ................................ ................................ ................................ ............... 59 Figure 4.5 : HWHM ( ) as a function of . Red line is fitted to Fickian model and blue line is fitted to Singwi - Sjolander model. ................................ ................................ ................................ . 60 Figure 4.6 : Nuclear density may of Cu atom showing the diffusion pathway: (a) 1 2 d 24g path; (b) 12d 12d path; (c) 12e 24g path. Brown atom is Sb and yellow is S. Isosurface level of 0.12 Å - 3 . ................................ ................................ ................................ ................................ ................ 61 xiii Figure 4.7: Diffusion paths schematic of (a) 24g 12d and (b) 24g 12e e xt racted from the actual MD trajectory and migration energy barriers of (c) 24g 12d and (d) 24g 12e from NEB calculations. ................................ ................................ ................................ ......................... 62 Figure 5.1: Average structure after 1x1x1 PW NVE simulation sh owi ng the crystal structure of Cu 10 Zn 2 Sb 4 S 13 . ................................ ................................ ................................ .............................. 68 Figure 5.2 : Lattice parameters ( a , b , c ) of (a) 1x1x1 PW and (b) 1x1x1 AO cell as a function of time in NPT simulation. ................................ ................................ ................................ ................ 68 Figure 5.3 : Partial and total vibrational density of states (VDOS) of each atom group (a - e: Cu 12e , Cu 12d , Zn 12d , Sb, and S; f: total) for three simulations (1x 1x1 PW, 1x1x1 AO, and 2x2x2 AO). ................................ ................................ ................................ ................................ ... 69 Figure 5.4 : Heat capacity per atom for three simulations (1x1x1 PW, 1x1x1 AO, and 2x2x2 AO). Experimental data from Lu et al. 7 and Lara - Curzio et al. 33 are also shown for comparison. ................................ ................................ ................................ ................................ ....................... 71 Figure 5.5 : Longitudinal phonon dispersion for (a) Q=0.60 A - 1 , (c) Q=0.84 A - 1 and transverse phonon dispersion for (b) Q=0.60 A - 1 , (d) Q=0.84 A - 1 . ................................ ................................ 72 Figure 5.6 : k 2 - weighted Cu EXAFS spectrum obtained by experiment (black circle) and simulations (solid lines). ................................ ................................ ................................ ............... 73 Figure 5.7 : Longitudinal (a) and transverse (b) response of the momentum correlation . (c) Di spersion of longitudinal (L) and transverse (T) sound modes. Peak energy for each Q value (red dots: longitudinal; blue dots: transverse) and their linear regression have been shown. ...... 75 Figure 5. 8 : Longitudinal (a) and transverse (b) response of the coherent velocity correlation. Longitudinal acoustic (LA) and quasi - localized (QL) modes are marked in a. The dotted lines have the same energy between a and b. ................................ ................................ ........................ 76 Figure 6.1: Flowchart of generation of amorphous a - Li 3 PO 4 and a - LiPON. ................................ 80 Figure 6.2: Schematic of a - LiPON structure at 1200 K from NVT run. ................................ ..... 82 Figure 6.3: Bond order as a function of bond distance of P - N, P - O and Li - X pairs. .................. 83 Figure 6.4: ELF maps of (a) Li - N a , (b) Li - N d and (c) Li - O units (isosu rface level of 0.82). ...... 84 Figure 6.5: Statistical results of Li coordination number around N a , N d and O. ......................... 85 Figure 6.6 : An example of as a function of Q 2 in Li diffusing at 1200 K for a - LiPON. .......... 85 Figure 6.7: of P and O for c - Li 3 PO 4 and c - LiPON at 1400 K and 1500 K. ........................ 86 Figure 6.8: Self - diffusivity of P/O/N atomic groups from Li 3 PO 4 and LiPON as a function of temperature. ................................ ................................ ................................ ................................ .. 88 xiv Figure 6.9: Li self - diffusivity of various structures as a function of temp erature comparing to experimental results. ................................ ................................ ................................ ..................... 89 Figure 6.10: Residence time (a) and jump length (b) of Li diffusion. ................................ ......... 89 Figure 6.11: (a ) Time and (b) frequency domain for the transverse component of coherent charge current density function for three smallest Q at 1400 K for a - LiP ON. ................................ ......... 91 Figure 6.12: Calculat ed ionic conductivity of various structures as a function of temperature comparing to experimental data. ................................ ................................ ................................ ... 92 Figure 7.1: Initial supercell of LiPON/Li interface. ................................ ................................ ..... 95 Figure 7.2: Average structure over 10 to 80 ps during the NVT run. ................................ ........... 96 Figure 7.3: Snapshots of LiPON/Li interphase at (a) 10 ps, (b) 20 ps and (c) 80 ps. ................... 97 Figure 7.4: Variation of Li charge from DDEC charge calculation along the c axis. .................. 98 Figure 7.5: (a) A snapshot from MD at 20 ps with 3 s elected layers (thickness o f 6 Å) and two chosen Li atoms from metallic Li phase. (b) Projected density of s tates at selected layers and atoms. Fermi energy level is located at zero energy. ................................ ................................ .... 99 1 CHAPTER 1 : Introduction and motivatio n 1.1 Energy crisis Energy crisis is a global issue in the recent decades due to the decreasing supply and increasing demand of fossil fuels. According to the BP Stat istical Review of World Energy Report 1 , the overall energy consumption keeps increasing for the past years , as shown in Fig ure 1 .1 , due to the e xpress industrial development and robust global economy. It can be noticed from the recent 2018 data that the non - renewable energy sources such as natural gas , coal and oil, we re still playing an important role in the energy demand , wh ile the application of renew able energy resources contributed only 9.5% to the overall energy consumption, which wa s still limited and inadequate. Figure 1 . 1 : Global energy consumption from 1993 to 2018 wit h the un it of mtoe (million tonnes oil equivalent) 1 . 2 In addit ion, d riven by the higher energy demand, the energy - related CO 2 emission has continually increased over the past 20 years and reaches a historic high of 33.1 Gt in 2018 as revealed in Figure 1. 2 2 - 3 . The increase of CO 2 content in the atmosphere have been accused to be one of the main reasons of the global clima te change and thus lead to num bers of negative effects such as global warming 4 - 7 . Figure 1 . 2 : Global energy - r elated carbon dioxide emissions by source, 1990 - 2018 2 . The causes of the energy crisis are complex and can be summarized as overconsumption of natural resources, overpopulation, un developed infrastructures, unexplored renewable energy options , waste of e nergy, etc. Considering the enor mous amount of waste energy, recovery of those can reduce our energy consumption and improve the energy efficiency. According to statistical analysis from Forman et al. 8 , about 72% of the energy consumed from primary energy sources were wasted as rejected energy, mainly in the form of heat, while the rmoelectric can help to recover t hose waste heat by converting thermal energy to ele ctrical energy 9 . Another solution is moving towards renewab le resources. However, the renewab le resources, such as solar energy, 3 wind power and tide power, are intermittent, limiting their large - scale application. Thus, energy storage devices are necessary to store the excess energy. Lithium - ion battery is current ly one of the most successful and reliable energy storage devices, attracting many re search interests 10 . In this thesis, two promising energy materials, tetrahedrite Cu 12 Sb 4 S 13 as thermoelectric and Lithium phosphorous oxynitride ( LiPON ) as solid - state electrolyte for Li - ion battery , will be discussed. 1.2 Thermoelectric materials and te trahedrite Cu 12 Sb 4 S 13 One effective way to ease the energy demand and reduce our dependence of fossil fuel, is to utilize the huge amount of unused waste heat and convert them to electrical energy. Thermoelectric material s are promising due to their abilit y of direction conversion between thermal and electrical energy , which can be used to recover waste heat or solid - state coolers 11 - 13 . It has been reported that the re are 191 million vehicles in US, dissipating 66% of their consumed energy as wasted heat emission, producing about 36 TWh of waste process heat per year 1 4 . Considering the enor mous wasted heat s not only from automobiles, but also from industrial operations , human activities, etc., researchers are throwing their efforts to find high - efficient and reliable thermoelectri c materials . 1.2.1 Thermoelectric ef fect and performance of thermoelectric materials Thermoelectric effect is known as the direct conversion between temperature difference and electric voltage 15 . The phenomenon that thermoelectric devices can convert thermal energy into electrical power is k overed in 1821 by Thomas Seeback 16 propos ed in 1834 17 . Figure 1. 3 shows the schematic of thermoelectric device which is the application of thermoelectric effect . Generally thermoelectric device contains a lot of thermoelectric couples wired electrically in series and thermally in parallel, whi le each of the couple consists of an n - type 4 (negative thermalpower and electron carriers) and a p - type (positive thermalpo wer and hole carriers) semiconductor. In the power generator case, the temperature difference ( ) provides the voltage ( ) due to t he Seebeck effect ( ) , where is the Seebeck coefficient, and thus drives the electrical current and generates po wer output. In the cooling mode, the electrical current ( ) supplied by external power source leads to a heat flow ( ) from top to b ottom because of the Peltier effect ( ) and thus cool s the top panel 11 , 18 . Figure 1 . 3 : Schematic of thermoelectric device showing the direction of charge flow on both cooling and power generation 11 . 5 (a) (b) Figure 1 . 4 : (a) Optimizing through carrier concentration tuning. (b) Benefits of reducing . 11 The per formance of thermoelectric materials is quantified by a dimensionless figure of merit : , ( 1 ) where the electrical conductivity, the absolute temperature, and the thermal conductivity. Thermoelectric materials such as SnSe 19 and P bTe 20 could reach zT at or near the commercial threshold value of 2.5 at high temperatures, enabling the their application in the recovery of waste heat. Fig ure 1. 4 i llustrates the relationshi p bet ween and these parameters with carrier concentration tuning. Figure 1. 4 (a) shows that t he peak of locates between 10 19 to 10 20 carrier/cm 3 , representing hea vily doped semiconductor. It can be noted that Seebeck coefficient ( in the refer r ed figure ) and electrical conductivity are inversely related . Moreover, low thermal 6 conductivity is necessary for thermoelectrics from equation . Thermal conductivity has two components, electronic thermal conductivity and lattice thermal conduct ivity . However, and are correlated by the Wiedemann - Franz law : , ( 2 ) where is the Lorenz factor with 2.4 x 10 - 8 J 2 K - 2 C - 2 for free electron . This equation suggest s that simply inc reasing comes with the increase of and usually decrease of thermal power 9 . In terms of , Figure 4 (b) shows the benefits of having a low , where point (1) is the optimized with = 0.8 W m - 1 K - 1 and that is a function of carrier concentration for a model materials Bi 2 Te 3 . If we reduce to 0.2 W m - 1 K - 1 , will increase to point (2) while keeping unchanged. Furthermore, reducing by decreas ing carrier concentration can reoptimized the value of and reach the maximum at point (3) 11 . Considering all the relations and conflicts among these parameters , - glass electron - Slack 21 , suggest ing that the ideal thermoelectric s require a ra ther unusual material: a low lattice thermal conductivity as glass, and a high electrical conductivity as crystal. 1.2.2 Application s of thermoelectric materials Advantages of using solid - state thermoelectric materials include compactness, quietness (no moving parts inside), and localized heating and cooling 18 . It has been mentioned above that thermoelectric materials have power generation mode and cooling mode , accompanied with considerable ef fic iency, thus the application s of thermoelectric materials can be everywhere in either industry or personal lives . In terms of power generation, automobiles could be one of the main applications of thermoelectric materials since a lot of waste heat are ge ner ated during the exhaust of vehicles. By converting the waste heat into electrical energy using thermoelect rics, the efficiency of fuel can be 7 improved and thus reduce the emission of greenhouse gases 22 - 23 . Furthermore, i n th e field of power - generating plants, about 67% of its available energy are wasted during the electrical pro duction, while in manufacturing industries, about 33% of their energy are rejected to the atmosphere in the form of heat. These waste thermal ener gy a re equivalent to the energy generated by more than 1 billion barrels of oil 24 - 2 5 . In these aspects, thermoelectrics have great potentials to recover the waste heat and reproduce enormous amounts of e lectr ical power. In addition, thermoelectrics with high value can be applied in the case of converting solar thermal energy and replacing the currently used dynamic converter 26 . Cooling applications of thermoelectrics are also emerging and can be widely used in different areas 11 . For instance , thermoelectrics cooling devi ces can be commercialized as portable refrigerator, beverage can cooler and picnic basket in the civil market 27 - 29 . Moreover, electronic devices like PC processors o r cell p hones generate a large amo unt of heat while stable performance of these microelectronic s require s a relatively low temperature. In these case s , thermoelectric cooler attract s great attention due to their small size and quietness 30 . Besides, there has been great inte rest s in automobile cooling applica tions, medical applications, air - conditioning applications and industrial temperature control applications of thermoelectric coolers 31 - 33 . 1.2.3 Tetrahedrite Cu 12 Sb 4 S 13 obvious that low thermal conductivity is critical for thermoelectric materials performance and thus offer the guideline of research and selection of optimal materi als. Tetrahedrites, the copper antimony sulfosalt mineral, typified by Cu 12 - x M x Sb 4 S 13 , where M is a transition metal element such as Ni, Zn, Fe or Mn, ha ve the potential for thermoelectric application due to their e arth - abund ance, environmental friendliness, favorable electrical properties, and most importantly 8 intrinsic low lattice thermal conductivity (less than 1 W m - 1 K - 1 ) in wide temperature range 34 - 36 . Lu et al reported Zn - doped tetrahedrite Cu 12 - x Zn x Sb 4 S 13 reach low (< 0.5 W m - 1 K - 1 ) close to spacing , and also obtain a high value ~ 1 at 720 K 35 . Figure 1 . 5 : Crystal structure of tetrahedrite Cu 12 Sb 4 S 13 with two distinct Cu sites. Crystal structure of tetrahedrite Cu 12 Sb 4 S 13 , a body - centered cubic structure with I 3 symmetry, ha s been revealed in Figure 1. 5 by the software V E STA 37 . Two distinct Cu sites , identified as Cu 12d and Cu 12e , two S sites, S 24g and S 2a , and one Sb site, Sb 8c , are the components in the system . Cu 12d atoms are surrounded by four S 24g atoms wi th tetrahedra coordination , while Cu 12e atoms coordinated by t wo S 24g atoms and one S 2a atom and form a triangular planar coordination , and also bonded to two Sb atoms and form a Sb[CuS 3 ]Sb trigonal bipyramid . All the transition metal will take the Cu 12d s it es. Each S 24g is coordinated to two C u 12d atoms, one Cu 12e atom and one Sb 8c atom, whereas each S 2a is bonded to six Cu 12e 9 atoms in octahedral coordination. For the Sb site, each Sb 8c is coordinated with three S 24g atoms and three Cu 12e atoms. The origi n of low lattice thermal conductivity is strongly related to the structure of tetrahedrite . The temperature dependent single crystal XRD results performed by Pfitzner et al. 38 indicated a large anharmonic displacement of Cu 12e atom perpendicular to the plane of coordinating S atoms. Detailed study on the bonding environment, C u 12 e out - of - plane movement and Sb lone pair was carried out by Lai et al. 39 by first - principle simulation and synchrotron diffraction experiments, and linked the anharmonic rattling to the local bonding asymmetry and thus leaded to low thermal conductivity. Electronic behaviors of tetrahedrite is attracting many researchers. Lu et al. 35 calculated the band structure and density of states of Cu 12 Sb 4 S 13 and C u 1 0 Zn 2 Sb 4 S 13 using density - functional theory (DFT) , as shown in Figure 1.6. Results suggested that Cu 1 2 Sb 4 S 13 is metal - like structure, with Fermi energy level located near a sharp perk at top of valence band, while a semiconducting bandgap of ~ 1.1 eV can b e identified. The valence band s are primarily formed by hybridization of S 3p and Cu 3d orbitals. When Zn was doped into the structure as Cu 1 0 Zn 2 Sb 4 S 13 , the Fermi energy level was pushed towards the gap . Hence, the Zn - doped tetrahedrite behaved as a true s emiconductor , as shown in Figure 1.6 (b). Followed by this work, many researchers applied similar method on doping or substituting tetrahedrite in order to adjust the electronic properties 40 - 44 . 10 (a) ( b ) Figure 1 . 6 : Electronic band structure and density of states of (a) Cu 12 Sb 4 S 13 and (b) Cu 10 Zn 2 Sb 4 S 13 35 . Be s ides the Zn dopants, doping effects for various transition metal such as Mn, Fe, Co and Ni haven been studied over these years 45 - 46 . Appropriate dopants can not only adjust the electronic properties but also enhance the figure of merit, zT . One of the highest zT value ( zT =1 at 720 K) 11 was achieved by Lu et al. 34 by co - doping Ni and Zn into tetrahedrite. Recently , substitution on S b site by other elements, such as Te and Bi, has also been studied , with significant improvement to the performance of thermoelectrics 47 - 49 . Althoug h there are a lot of researches in the field of tetrahedrites including enhancement o f thermoelectric performance 34 , 36 , 50 , electronic p roperties 40 - 44 , different synthesis methods 51 - 53 , etc . , there are still plent y of un known s and uncertaint ies remained such as the dynamic properties, th e Cu - rich tetrahedrite system, and the unique low - temperature behavior s , which will be discussed in detail in Chapter 2, 3, 4 and 5 in this thesis . 1.3 Lithium - ion battery and elect rolyte material LiPON 1.3.1 Introduction of lithium - ion battery Recalled the energy crisis we discussed above, besides the help from thermoelectrics, renewable energy sources such as solar energy 54 - 56 , wind power 57 - 59 and bioenergy 60 - 62 are making great contributions. However, one issue which cannot be igno red is the mismatch between excess energy ge neration a nd energy consumption . Therefore, reliable energy storage devices are emerging to relieve the energy demand and improve the efficiency of consuming energy resources . Lithium - ion batteries are known to b e reliable and successful electrochemical en ergy stora ge devices and can found everywhere in our daily lives, including laptops, smartphones and electrical vehicles. The general working principles of Li - ion batteries are shown in Figure 1. 7 using Li x C 6 /Li 1 - x CoO 2 cell as an example. A Li - ion battery is compos ed of three essential parts, cathode (positive electrode), anode (negative electrode) and electrolyte (separator). During the charging process, lithium ions shuttle from the cathode to the anode thro ugh the electrolyte while the electrons move 12 fr om anode t o cathode through external circuit simultaneously. And the reverse process occur s when discharging. Figure 1 . 7 : Schematic of working principles of Li x C 6 /Li 1 - x CoO 2 Li - ion battery 10 . A suggested road map for development of electrode materials for Li - ion bat teries is showed in Figure 1. 8 . In terms of cathode materials, Li Co O 2 was firstly considered and applied in electronic devices. However, the disadvantages such as high cost, toxicity, safety issue and relatively low capacity (~135 m Ah g - 1 ) limit ed its larg e - scale application 63 . Later LiMn 1.5 Ni 0.5 O 4 was developed and proved to have high working voltage of 4.7 V (compared to 4.1 V for LiCoO 2 ). LiMPO 4 (M=Co , Fe, Mn, Ni or combinations of those) compounds are also becoming attractive due to their high capacity, e.g., 197 m Ah g - 1 for Li 3 V 2 PO 4 64 . Whereas, there are still challenges 13 like rate capability and cycling performa nce f or those to b e commercially used. Besides, Li - excess layered oxides such as (Li 2 MnO 3 ) x (LiMO 2 ) 1 - x (M=Ni, Co, Mn) offered higher theoretical capacity (~250 m Ah g - 1 ) with a working voltage of ~4.0 V 65 . It is notable fro m Figure 7 that cathode materials like Li - O 2 , Li - S/C or Li 2 S - Si could be promising cathodes since they could offer even higher capacity. However, reversibility, compatibility and cycling performance are still big issues and need to be improved in the futur e 66 . Figure 1 . 8 : Current electrode materials road map 67 . On the anode side, carbon in the form of graphite is the most widely used anode material in Li - ion cells due to it s capability of intercalation and deintercalation of Li ions into and fr om its layered s tructure 68 - 69 . Although high energy density can be achieved by combining graphite and suitable cathodes, the low react potentials will easily lead to the formation of a passivating solid electrolyte inter phase (SEI) 63 on the graphite side and result in poor cycle performance of the cell. 14 Furthermore, volumetric expansion or shrinkage of graphite occur red during the charging and discharging processes will p ossibly reduce the inter facial contact and lead to capacity fading. extremely small volumetric change ( ~ 0.2%) during Li ion intercalation and deintercalation 70 . As a result, LTO has the characteristics of good cycle life , high rate capability and excellent safety features. While graphite can only accommodate one Li ato m per graphene unit and t hus having theoretical electrochemical capacity of 372 mAh g - 1 , other metals and metalloids such as Sn and Si, can accommodate more than four Li atoms per unit and thereby reaching much high capacities of 960 and 4009 mAh g - 1 , resp ectively 71 . Figure 1 . 9 : Tem perature dependence of ionic conductivity of various electrolyte materials 72 . Turning to electrolyte materials, they can be organized into several groups: inorganic solid - state electrolytes, organic l iq uid electrolytes, pol ymer electrolytes, ionic liquid and gel electrolytes. 15 Because electrolytes are primarily used to conduct ions during charging and discharging processes, ionic conductivity is one of the most critical criteria when selecting suitable ca ndidates. The tempera ture dependence of ionic conductivity of commonly used Li electrolytes has been summarized in Figure 1. 9 . Among these various kinds of electrolyte materials, liquid electrolytes such as LiBF 4 , LiPF 6 dissolved in organic solvents like e thylene carbonate (EC ) and dimethyl carbonate (DMC), have the disadvantages of poor electrochemical stability, limited operation temperature range, leakage, and most importantly the safety concerns (flammable) 73 . The polymer electrolytes, although having characteristics of great safety, easily preparat ion and flexible shape, there are many probl ems existed including instability of electrolyte/electrode interface, narrow operation temperature and limited mechanical properties 74 - 76 . Thus, in order to select proper electrolyte mat erials for large - scale application without safety issues, inorganic solid - state electrolytes might be the solution. Firstly, i norganic solid - state electrolytes display comparable i onic conductivity with liquid electrolytes , e.g., gar net - type oxides 77 exhibit high Li - ion conductivity of 10 - 3 S cm - 1 10 GeP 2 S 2 72 offers superior Li ionic conductivity of 1.2 x 10 - 2 S cm - 1 at room temperature. Secondly, solid electrolytes generally perform wider electrochemical window than liquid electrolyte s , which enable the cell to be operated under high cut - off voltage. In ad dition, good thermal stability of solids widens the operation temperature of batteries. Last but not the least, great safety operating environment can be achieved. 1.3.2 Electrolyte m a terial LiPON Lithium phosphorous oxynitride ( LiPON ) with general formula Li x PO y N z , where x=2y+3z - 5, are widely used as thin - film solid - state electrolytes in Li - ion battery. They were pioneered at Oak Ridge National Laboratory by sputtering Li 3 PO 4 target i n N 2 gas, while the ionic conductivity increased fr om 7 x 10 - 8 S cm - 1 to 2. 2 x 10 - 6 S cm - 1 at room temperature 78 . Figure 1. 10 16 reveals the quaternary diagram of stoichiometries of LiPON with four starting compounds ( LiO 1/2 , LiN 1/3 , PO 5/2 and PN 5/3 ). M ost of the reported stable or meta - stable phases among t he LiPON family are included. Besides the stoichiometric relationships, various structural patterns can be observed among the LiPON fa mily. The tetrahedral PO 4 - x N x can be identified for all LiPON structures, but with different connecting options 79 - 82 . It can be [ P O 4 ] - 3 in Li 3 PO 4 , [PN 4 ] - 7 in Li 7 PN 4 , or phosphate dimers [P 2 O 7 ] - 4 in Li 4 P 2 O 7 compound. In addition, the long chains pattern by phosphate groups sharing tetra hed ral corner atoms can be found in LiPO 3 83 . Figure 1 . 10 : Composition diagram with starting materials of LiO 1/2 , LiN 1/3 , PO 5/2 and PN 5/3 . Natu ral and synthetic crystalline materials are labeled as dark (blue) circles, thin - film compositions as light (turquoise) circles, and examples of stable and meta - stables nitride phosphate materials as black square 83 . Another interesting aspect in the LiPON family is the local environment of doping N atoms. Singly coordinated N can be found in Li 7 PN 4 , while doubly coordinated N was observed in 17 LiPN 2 84 - 85 . Moreover, N can b e both doubly and triply coordinated in structures like - P 3 N 5 86 . However, th e actual local environment of doping N and its coordination in LiPON remains inconclusive and will be discussed in the following chapter. LiPON is one of the most frequently used solid electrolytes in thin - film batteries. It had been reported that thin fil m batteries with LiCoO 2 as cathode, Li metal as anode and LiPON as e lectrolyte, are quite stable during the charge - discharge reactions (20,000 cycles with 0.001% capacity loss) 87 - 88 . Although t he ionic conductivity of LiPON is relatively low comparing to other solid - state electrolytes in Figure 1.8, it is the only demonstrated solid - state electrolyte which is quite stable in direct contact with Li metal at potentials from 0 - 5 V , as shown in F i gure 1. 11, which enable s the potential of doubling the energy density of current commercial batteries 89 - 91 . In addition, LiPON can be used as additional stabl e pro tective SEI on top of liquid or polymer electrolyte to prevent reaction between electrodes and electrolyte 92 . How ever, the structure of LiPON , the effects of N doping, the interface reaction at LiPON /Li and the origin of its good electrochemical stability remains inconclusive , which we will discuss more in Chapter 6 and 7 . Figure 1 . 11 : Current - potential curve of a LiPON / P t half - cell with respect to Li/Li + reference 91 . 18 1.4 Modeling techniques on structural and dynamic al propertie s of energy materials Considering the current experimental techniques have limits on studying specific research problems, computational modeling methods , which construct mathematical models to numerically study the behaviors of complex systems by simulatio n of computers, are proven to be quite powerful , efficient and reliable in many subjects including materials science, physics, chemistry, biology, etc . Figure 1. 1 2 shows some popular computational techniques for different length and tim e scales. In this ch apter we will discuss d ensity - functional theory (DFT) and density - functional tight - binding (DFTB) in details since they are the computational techniques we will include in the following chapters for tetrahedrites and LiPON systems. Figure 1 . 12 : Popular computational techniques for different length and time scales. 1.4.1 Introduction to DFT DFT, developed by Hohenberg, Kohn and Sham 93 - 94 , is a computational quantum mechanical modelling method to investigate the electronic structure of many - body syst e ms . DFT locates in the bottom left in Figure 1.9 indicating it is applied in the smallest length and time scale, 19 while it is the most accurate one as well. DFT reformulates the Schrödinger equation (SE), which describes the behavior of electron in a syste m , so that analytical solutions to SE can be t ractable for real materials. Kohn and Sham 94 defined the most complex electron interac tions in an exchange - correlation functional, whose e xplic it form remains unknown and need to be approximated. T hu s, the selection of exchange - correlation functional, in a way, determines the accuracy of DFT calculation 95 - 96 . The foundation of modeling material properties is based on the interactions between the electrons and nuclei, which follow the fundamental laws of quantum mechanics (QM). The basic of Q M i s the solution of SE, with the form of: , ( 3 ) where is wave function of the system, and is a Hamiltonian and an eigenvalue of the wave function, respectively. However, in practical cases, it is difficult to s olv e the exact i n many - electrons systems , thus, the Born - Oppenheimer (BO) approximation 97 has been made to solve SE by considering independent motions of electrons and nuclei. Nonempirical a b initio simulation with BO can calculate the system energy accurately, however, wit h ve ry expensive computational cost and hence limits the numbers of ato ms and the timescale. Later some empirical or semi - empirical methods are developed to research the material properties. Among them, DFT is the most accurate and least expensive, and hen ce t he most popular one. The fundamental of DFT based on two theorems by Hohenberg 93 and Kohn 94 : Hohenberg proposed that the ground - state wave functi on is a unique functional of the electron density function and th erefore making the SE tractable by reducing the degree of freedom in SE from 3N to 3 dimensions; Kohn presented that the defined energy functional can be minimized by a 20 ground - state el ect ron density . With these two theorems, the total energy within DFT can be expressed as: , ( 4 ) where is the ki net ic energy of a noninteracting system with density , is the external potential of electron and neutron interactions, the third term refers to the potentials of electron - electron interactions, and the last term corresponds to the electron exchange - correlation energy of an interacting system with den sity . Kohn and Sham 94 also showed that the charge density can be represented by: , ( 5 ) , ( 6 ) w here is the wave function of electron for SE , is the effective external potential, is the eigenvalues. As we mentioned above, the exact form of is still unknown. However, if varies sufficiently, the can be written as: , ( 7 ) whe re is defined as exchange and correlation energy per electron with density . The approximation of exchange - correlation energy is one of the most important steps during practical DFT application s , which has been developed by many researchers. Lo cal density approximation (LDA) 98 or local - spin - density approximation (LSDA) 99 are the most successful ones, which have be en pro ven to b e reliable for solids. While generalized gradient approximation (GGA) is more accurate by applying gradient corrections, it is getting more and more attractions . Perdew - Burke - Ernzerhof (PBE) parametrization 100 within GGA is now the most widely used 21 model, especially in materials science. Later on, hybrids method was proposed by replacing a fraction of GGA exchange wi th Hartree - Fock (HF) exchange, fo r examp le, B3LYP method , which is more popular among chemists 101 . Figure 1.1 3 shows the number of recent DFT citations using these two method s, PBE and B3LYP. Besides, the Dark legend in the figure refers to those papers using either one without citing the orig inal papers, and Other means all other DFT - related papers. Figure 1 . 13 : The number of DFT citations from 1995 to 2013 102 . Molecular dynamics (MD) is wide ly used combin ing with DFT method, which describes the motion of collected atoms by their p ositions and momenta , and estimate the macroscopic physical properties by mean s of stat istical mechanics 103 . Time evolution of the system can be , ( 8 ) where is the force acting on atom at time , the atomic mass, the position , and the interaction potential. A practical MD simulation r equires an initial unit cell, periodic boundary condition, proper selection of interaction potentials, and stabilized controlment of temperature and 22 pressure. MD simulation can be accomplished in various ensembles, suc h as grand canonical ( VT), canonical (NVT), microcanonical (NVE) and isothermal - isobaric (NPT). The temperature and pressure can be adjusted by pre - defined thermostat (e.g., Langevin, Berendsen, Nose - Hoover) and barostat (e.g. Parrinello - Rahman, Andersen, Berendsen), respectively. The general inputs of DFT include the coordinates of atoms, boundary condition of unit cell, the choice of exchange - correlation functional, the parameters and algorithms for numerical and iterative convergence, and the way to tre at electron system (typically the use of pseudopotentials). The outputs of DFT are useful quantities such as the total energy of system, cohesive energies, energy barriers, charge density of electrons, electronic band structure, etc 104 - 105 . DFT is proven to be po werful in predicting properties of energy materials and des i gn of new materials in various aspects, such as thermoelectrics 39 , Li - ion batteries 106 - 107 , hydrogen production and storage 108 , superconductors 109 , photovoltaics 110 , etc 111 . Therefore, in this thesis we apply the reliable and accurate DFT method to study the structural, electronic and dynamics pro perties of tetrahedrite the rmoelectrics in Chapter 2, 3, 4 and 5. 1.4.2 Introduction to DFTB However, the limitation s of DFT still exist. For example, in a highly correlated electron systems, most functionals will fail. In addition, considering current com putatio nal expense, the sim ulated cell size cannot be too large (typically 1,000 atoms or less) and the time period should be short, th erefore, long - time dynamics cannot be extracted. Other limitations include the difficulty in modelling va n del Waals inte raction s and non - ground sta te properties 112 . In recent years, DFTB method was derived from DFT by neglection, approximation and parametrization of inte grals 113 - 115 , with the advantage of 2 - 3 orders of magnitude fas ter in computational speed th erefore allowing the treatment of larger unit cell and longer simula tion time . 23 The starting point of DFTB based on the second order expansion of the total energy expression of DFT 113 . The original form of total energy within DFT can be found in Equation (4). In a system of M electron s, N nucle i, at position , it may be re written as: , ( 9 ) w here the first term refers to over occupied Kohn - Sham eigenstates , the second is the exchange - correlation contribution and the third one is ion - ion core repulsion, . Then this equation can be re formulated by substituting the charge density by a sup erposition of neutral atomic density , , and expandi ng at the reference density to the second order in the density fluctuations, as the following: , ( 10 ) where represents the . From here, several app roximations have been made to construct the DFTB total energy expression: 1. In the second - order term, the charge density fluctuation is represented by , while can be appro ximated by atomic charge fluctuations , where is determined by Muliken population analysis and is the number of electrons for neutral atom . A function , depending on the Hubbard parameter ( is the second derivative of the total energy of a single atom with respect to the occupation number of the highest occupied atomic orbital) , is introduced to represent the integral over term and second derivative of . Therefore, the second - o rd er term can be expressed as: 24 . ( 11 ) 2. Kohn - Sham eigenstates are expanded in a minimal basic set , which consist s a linear combination of atomic orbitals 116 - 117 , in the form of: . ( 12 ) Here is determined by solving a modified Kohn - Sham equation: ( 13 ) where is the kinetic energy operator and the effective potential of a f ree neutral pseudoatom. Originally , the compression radius, was treated as a variational parameter 116 but in practice it is ch osen to be 2 times the covalent radius 113 . Here has a default value of 2. When two - center approximation has been adjusted for non - diagonal elements and eigenvalues of the free atom serve as diagonal elements, t he Hamiltonian matrix elements yielded as: . ( 14 ) The Hamilton can be precalculated and tabulated with the overlap matrix for interatomic distances and leads to a new expression as: , ( 15 ) where . The solution of the eigenvalue pr oblem for dominates the computational cost of DFTB calculation and hence provide the efficiency of this method. 25 3. For the r emaining terms of Eq. (10) , they can be summarized and defined as a short - range repulsive energy as: . ( 16 ) is then approximated as a sum of two body potentials 118 as: , ( 17 ) where refers to atom - type specific pair potentials , and is constructed as the difference between the total energy calculated in DFT and the corresponding electronic energy obtained by DFTB for the proper systems 114 , 119 as: . ( 18 ) Thus, combining the approximations and substitutions discussed above, DFTB total energy formalism shows as: , ( 19 ) where the first term contains the DFTB matrix elements and the third term describes the repulsive potential. These two terms together refer to standard (non - sel f - consistent) DFTB method 117 , 120 , and the second term corresponds to the self - consistent charge (SCC) DFTB by approximation of the second order term of DFT Taylor series expansion. DFTB, comparing to DFT, is efficient and s uitable to study larger cell (more sampling) and long - time scale simulation, hence being applied to study the structure and dynamic properties such as diffusivity and ionic conductivity of LiPON in Chapter 6. We al so use this method to investigate LiPON / Li metal interphase in Chapter 7 since a much larger cell is required. 26 1.4.3 Other computational techniques There are many other available and successful modeling techniques in energy materials simulation , for example, i nteratomic p otentials. I nteratomic po tentials use mathematical functions to calculate the potential energy of a system given atomic position in space 121 , which are popular due to its capability treating a large system (thousands of atoms). General interatomic potentials expression can be written as a series expansion of functional terms that depend on the posit ion of atoms, while the total energy of the system becomes 122 : , ( 20 ) where , and is the one - body term, two - body term and three - body term, respectively, N the number of atoms, the position of atom i . The forces defined between atoms then can be obtained by taking three - dimensional derivative of with respect to the position of atom i as: . ( 21 ) The interatomic potentials have various kinds when dealing with different physical objects, including pair potentials (for example, Leonard - Jones potential 123 - 124 ), many - body potentials (embedded atom model 125 , the Stillinger - Weber potentials 126 ) and short - range repulsive potentials. Although interatomic potentials metho d locat es in the top right in Figure 1.9, it is still under the category of molecular scale methods, as well as DFT and DFTB. In addition to molecular scale methods, there are various modeling techniques 103 dealing with the simulation at microscale such as Brownian dynamics (BD) 127 - 128 , dissipative particle dynamics (DPD) 129 , lattice Boltzmann (LB) 130 , tim e - dependent Ginzburg - Landau (TDGL) method 131 and dynamic DFT method 132 . Furthermore, there are mesoscale and macroscale methods , suc h as mic romechanics 133 and fini te element analysis (FEA) 134 - 135 , to predict and analyze the macroscopic properties of materia ls . To 27 sol ve a resear ch problem, it is helpful to combine two or even more modeling and computational methods discussed above to achieve a balance of efficiency and accuracy. 1.5 Dissertation outline This dissertation is organized in 8 chapters as the foll owing layo ut: In the present chapter, Chapter 1, the fundamental scientific backgrounds are included , such as the current status of energy crisis, thermoelectric effect , application and performance of thermoelectric materials and introduction of tetrahedri te thermoe lectric Cu 1 2 Sb 4 S 13 . Introduction to Li - ion batteries is also presented from the material science perspective, especially the electrolyte material LIiON . In Chapter 2, n eutron d iffraction and c omputational s tudy on pristine t etrahedrite Cu 12 Sb 4 S 13 is propos ed , focusing on it s l ow - t emperature structure and p roperties . A metal - semiconductor transition (MST) occurred at 85 K was observed by significantly increasement of electronic resistivity upon cooling, accompanied with abrupt changes of magnetic p roperties and thermoelectric properties as Seebeck coefficient and thermal conductivity 136 - 137 . However, the discussion that whether there is a cubic to tetragonal phase transition at the MST temperature remains inconclusive 138 . Inspired by these works, we perfo rm neutron d iffraction measurement and computational simulation on t etrahedrite Cu 12 Sb 4 S 13 . No cubic to tetragon al phase transition at MST temperature can be identified, however, negative thermal expansion below 50 K is observed, according to both experime ntal and com putational results. Sb - Cu 12e - Sb interaction is examined at different temperatures. Thermoelectric re lated properties such as DOS, electrical resistivity, Seebeck coefficient and electrical thermal conductivity are calculated comparing to experi mental data. 28 In Chapter 3, the incoherent and coherent dynamics study of tetrahedrite Cu 10. 5 NiZn 0.5 Sb 4 S 13 by DFT - based MD simulations have been proposed. Recent inelastic neutron scattering (INS) experiment indicated the anomalous behavior of tetrahedrite on Fourier transform of velocity autocorrelation function from simulation trajectory is presented comparing to INS data , where good consistency can be obs erved. Parti al VDOS for Cu 12e atoms at different temperatures shows that the low energy peak shifts to hi gher energy as temperature increases. The interaction between Cu 12e and Sb at different temperatures is also included, where we discover that the reduc tion of anha rmonic rattling of Cu between two Sb atoms upon cooling is responsible for the phonon softeni ng. Coherent atomic dynamics are also examined by computing the longitudinal and transverse dynamics structure factor based on the momentum correlation function. D ebye temperature, Debye heat capacity and sound velocity are calculated comparing to experime ntal data. Thermal conductivity is also extracted by making the approximation that the mean free path equals to the average interatomic spacing. In Cha pter 4, the transport properties study of Cu among Cu - rich tetrahedrite is presented using the Cu 14 Sb 4 S 13 as a model structure. Previous research 139 indicated that Cu becomes mobile above 393 K and are like ly to be superionic conductors in Cu - rich tetrahedrite , while the excess Cu atoms were like ly to be at Cu 24 g sites, however, the insight of distribution and movement of Cu atoms remains unknown. Incoherent density correlation is calculated, fr om which the vibrating Cu group and diffusing Cu group can be identified. Cu self - diffusivity, jump leng th and residence time are investigated. The nuclear density map is proposed to visualize the movement of Cu atom between different sites, where several diffusion pat hway can be observed. 29 ( NEB ) method is applied to estimate the migration energy barrier of different migration paths, where we select two actual paths ( 24g 12d and 24g 12e ) from trajectory as examples. In Chapter 5, we study the ef fect of cell size and basis set on results from DFT - based first - principles MD simula tion of a model thermoelectric material: Cu 10 Zn 2 Sb 4 S 13 tetrahedrite. Two cell sizes, 1x1x1 with 58 atoms and 2x2x2 with 464 atoms are selected, with two popul ar basis sets , plane waves (PW) and atomic orbitals (AO) . Structural and dynamic properties, such as lattice parameters, partial and total vibrational density of states, phonon dispersion, heat capacity, and EXAFS spectra , are extracted from each of the si mulation. Goo d agreement can be achieved by comparing with each other, except some minor differences, e.g., the partial VDOS, which we ascribe to different electronic structures and atomic forces from different basis sets. The larger 2x2x2 cell allow us to access the r egion closer to the hydrodynamic limit of atomic dynamics to study materials properties. In Chapter 6, structure and ionic conduction study on Li 3 PO 4 and LiPON with DFTB method is presented. N deposition into Li 3 PO 4 can increase the ionic cond uctivity from 7 x 10 - 8 S/cm without any N content to 3.3 x 10 - 6 S/cm according to previous study 91 , 140 . However, the origin of the N doping effect, e.g., the N coordination with P atoms, remains uncertain. Moreover, the study on the amorphous s tructure of LiPON is limited. Here in this chapter we apply the DFTB method to generate four model structures, crystalline Li 3 PO 4 ( c - Li 3 PO 4 ) , amorphous Li 3 PO 4 ( a - Li 3 PO 4 ) , crystalline LiPON (c - LiPON ) and amorphous LiPON ( a - LiPON ) , to investig ate their struc tural properties. Lattice parameters of crystalline phases are in good agreement with experimental data. Apical N (N a ) and doubly - coordinated N (N d ) can be identified in a - LiPON , whereas no triply - coordinated N (N t ) exists. Bond order (BO) c alculation is i mplemented to investigated different bonding environment around P and Li atoms. E lectron localized function 30 (ELF) maps and statistical analysis are performed to study the Li environment as well. In addition, we exa mine self - diffusivity of di fferent atomic groups . We discover that the addition of N in c - LIPON helps the diffusion of P and O in the crystalline phase by investigating the incoherent density correlation plots. Li self - diffusivity extracting by DFTB shows good consistency with exper imental results . Ionic conductivities for all four structures are also included by extrapolation of a Fourier transform form of transverse coherent charge current density correlation function from MD trajectories. In Chapter 7, p reliminary study on LiPON /L i interphase is proposed. A LiPON /Li interphase supercell is generated. No decomposition of LiPON , and no Li ions exchange between the bulk and metallic Li layers can be observed. A 6 Å interface layer is identified by investigating the Li charge distribut ion, beyond whi ch the metallic Li maintain as bcc structure. Projected DOS of selected layers and atoms are also revealed in this chapter, suggesting the metallic behavior of the interf ace. In Chapter 8, future works focusing on thermoelectric materials te trahedrite, LiP ON /Li interface and promising cathode materials are included, as well as the conclusions of the dissertation. 31 CHAPTER 2 : Pristine tetrahedrite Cu 12 Sb 4 S 13 : n eutron d iffraction and c omputational s tudy on l ow - t emperature p roperties 2.1 Introduction The wide applic ations of thermoelectric materials have been discussed above. Some important applications such as cryogen - free cryo genic solid - state cooling , required thermoelectric candidates with high efficiencies at low temperatures. However, study on th ermoelectric pr operties at low temperature region (below 300 K) is still limited. It has been reported that a metal - semiconductor transition (MST) at around 85 K will occur in Cu 12 Sb 4 S 13 46 , indicated by significant increasement of electronic resistivity, from about 2 x 10 - 5 10 - 4 at 50 K upo n cooling. Therm oelectric properties such as Seebeck coefficient and thermal conductivity, are also reported to have abrupt change at the same temperature 136 . However, it is still controversial that whether a structural phase transition accompanies the MST. For example, May et al 138 . reported a cubic to tetragonal phase transition below 88 K, suggested by cubic Bragg peak (611) splitting from powder x - ray diffraction. They also found out that Zn - doping would suppress the structural distort ion by comparing the tempe rature dependence of lattice parameters and phonon DOS from inelastic neutron scattering ( INS) measurements for Cu 12 Sb 4 S 13 and Cu 10 Zn 2 Sb 4 S 13 . Tanaka et al. compared synchrotron x - ray diffraction data of Cu 12 Sb 4 S 13 at 100 K and 75 K, and observed that (400) peak at 100 K spl it into double peaks at 75 K and concluded that the cubic to tetragonal phase transition was associated with MST 137 . On the other hand, Nasonova et al. did not observe structural transformation, nor any sign of reflection peaks splitting, by performing high - resolution power X - ray diffraction 136 . Another characteristic phenomenon that occurr ed in tetrahedrite, is a brupt drop of magnet ic susceptibility upon cooling, which also ascribed to the distribution of Cu 12e 137 , 141 . 32 In this work, we present crystal structu re study of stoichiome tric Cu 12 Sb 4 S 13 by neutron diffraction measurement in the temperature range of 20 300 K, comparing to experimental synchrontron x - ray data to have better understanding of the proposed structural phase transition. Density functional theory - based first - pri nciples molecular dynamics simulation is also applied to investigate the origin of the unique low - temperature features of Cu 12 Sb 4 S 13 and calculate its thermoelectric properties at low temperature range. 2.2 Computational and experimen tal details Vienna Ab initio Simulation Package ( VASP ) 142 - 145 code employing the Projector Augmented - Wave (PAW) method was used for molecular dynamics (MD) and Bo ltzWann 146 - 147 module of Wannier90 148 code. V alence electron configu rations for Cu, Sb, S ato ms are 4s 1 3d 10 , 5s 2 5p 3 , and 3s 2 3 p 4 , respectively. The plane wave energy cutoff was 450 eV. The generalized gradient approximation (GGA) with Perdew - Burke - Ernzerhof (PBE) parametrization 100 wa s applied as exchange - correlation (XC) functional and 2 x 2 x 2 k - mesh was used. Simulation cell was a 1 x 1 x 1 scale of unit cell with 58 atoms. For MD simulations, first a constant number of part icles, pressure and temper ature (NPT) ensemble was perform ed at 20, 50, 87.5 and 300 K for 3 ps to obtain lattice parameters. Then a constant number of particles, volume and energy (NVE) ensemble was implemented for 5 ps for equilibration and 20 ps for pro duction at same temperatur es as NPT. Time step was 1 fs fo r all the MD runs. Thermoelectric properties were obtained by BoltzWann module of the Wannier90 package using results from DFT - based MD as input 146 - 147 . Relaxation time of BoltzWann calculation is set to be 10 fs. Tetrahedrite (Cu 12 Sb 4 S 13 ) powder was synthesized by t he modified polyol process pr oposed by Weller et al 149 . To ensure uniform particle size, the powder was loaded into a stainless - 33 steel ball mill jar and milled for 10 min via a SPEX 8000D vibratory ball mill. Phase - puri t y of tetrahedrite powder was confirmed by X - ray diffraction (XRD) via a Rigaku Miniflex benchtop dispersive X - ray spectroscopy (EDS) was performed with a Hitachi TM - 3000 Tabletop Microscope, a Bruker XFlash MIN SVE detector, and an accelerating voltage of 15 kV. Data for atomic ratios were averaged over 4 sites and normalized to the measured S content. The measured composition according to EDS was Cu 13.4±0.3 Sb 4.1±0.1 S 13.0± 0.3 . Neutron diffraction was conducted at the BL - 11A POWGEN beamline, Spallation Neutron Source (SNS), Oak Ridge N ational Laboratory, using our synthesized Cu 12 Sb 4 S 13 sample, at 20 K, 50 K, 87.5 K, and 300 K. The software Jana2006 was used for structure refinement 150 . 2.3 Results and Discussion 2. 3.1 Low temperature structure of Cu 12 Sb 4 S 13 Neutron diffraction was conducted at four different temperatures as shown in Figure 2 .1 . Earlier reports 137 - 138 proposed that 004 peak was split into two peaks under MST temperature and thus Cu 12 Sb 4 S 13 underwent phase distortion f ro m cubic to tetragonal symmetry. In our case, no peak splitting of 004 reflection can be observed below the MST temperature, indicating that no structural phase transition occurred, and Cu 12 Sb 4 S 13 stayed as cubic phase, which i s in good agreement with the s ynchrontron x - ray data from Nasonova et al 136 . 34 Figure 2 . 1 : Neutron diffraction data for Cu 12 Sb 4 S 13 at 20, 50, 87.5 and 300 K, synchrontron x - ray data by Nasonova el al. are shown for comparison. No pea k s plitting can be identified below MST temperature. High quality of Rietveld refinement can be obtained for all temperatures using Cu 12 Sb 4 S 13 as a model as Figure 2 .2 and 2 .3 . Occupancy of Cu 12e and Cu 12d were relax ed, while occupancy of S and Sb were fix ed. Anisotropic atomic displacement parameters were used for all atom types. Atomic parameters at 20 K, 50 K, 87.5 K and 300 K are listed from Table 5. 1 to Table 5.4 . We also used a Cu - rich tetrahedrite model Cu 13.2 Sb 4 S 13 139 , for Rietveld refinement at 87.5 K for comparison, since this composition is close to our measure composition by EDS. 35 (a) (b) (c) (d) Figure 2 . 2 : Neutron diffraction data collected for tetrahedrite sample at 87.5 K using (a, b) Cu 12 Sb 4 S 13 and (c,d) Cu 13.2 Sb 4 S 13 as model structures. Block 1 contains 80% of the data and block 2 has t he rest 20%. 36 (a) (b) (c) (d) (e) (f) Figure 2 . 3 : Neutron diffraction data collected for tetrahedrite sample at (a, b) 20 K, (c, d) at 50 K and (e, f) at 300 K using Cu 12 Sb 4 S 13 structure as model. Block 1 contains 80% of the data and block 2 has the rest 20%. 37 Table 2 . 1 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 20 K. Table 2 . 2 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 50 K. Table 2 . 3 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 87.5 K. Table 2 . 4 : Atomic parameters in the crystal structure of Cu 12 Sb 4 S 13 at 300 K. Lattice parameters from DFT - MD NPT ensemble have been plotted in Figure 2 .4 , together with lattice parameters obtained from our neutron diffraction experiment and s ynchrontron x - ray from Nasonova 136 . MD simulation showed that Cu 12 Sb 4 S 13 maintained as cubic symmetry at the designed temperature range, which is consistent with our neutron experiment. MD data are about 0.4% higher than neutron and synchrontron data, indicating good relia b ility of our 38 simulation. For both MD and neutron lattice parameters, positive thermal expan sion from 50 300 K, and negative thermal expansion from 20 50 K can be observed. In comparison, synchrontron data have an upturn at 70 K upon cooling then reach e s local maximum at 40 K and decrease again. Previous report 136 suggested the negative thermal expansion is associate d with the position of Cu 12e atoms. Inspired by this, we analyzed MD simulation trajectory to explain this phenomenon. Figure 2 . 4 : Lattice parameters of Cu 12 Sb 4 S 13 obtained by DFT - MD NPT simulation and neutron diffraction refinement. Lattice parameters proposed by Nasonova using synchrontron x - ray are shown here for comparison. 2. 3.2 Sb - Cu 1 2e - Sb interaction To gain the atomic insight of the negative thermal expansi on and crystal structure of Cu 12 Sb 4 S 13 at low temperature, we inspected the distance between each Cu 12e and the two coordinated Sb atoms as the schematic in Figure 2 .5 (a). Accordin g to our previous research 39 , the rattler atoms (Cu 12e ) are regulated by the nearby lone pair electron of Sb and responsible for the 39 interested anharmonic movement, and lead to low lattice thermal conductivity of Cu 12 Sb 4 S 13 . Here we plotted the one - pair - potentials (O. P r.P) of Cu 12e - Sb by converting the distribution of dist ance between those two atoms at 20 K and 300 K. At 20 K, a double - well shape (highlighted red in Figure 2 .5 (b)) can be identified for the Cu 12e and Sb interaction energy, which indicated that the associated Cu 12e a tom is locked to one of the bonded Sb with average distance of 2.95 Å. 6 of total 12 Cu 12e atoms in our unit cell appeared to be locked with bond order of 0.2169 based on the Density Derived Electrostatic and Chemical scheme 151 calculation, while the other 6 12e are rattling between two Sb atoms, as s h own in Table 2 .5 quartic shape (highlighted blue in Figure 2 .5 (b) and (c)) can be observed and suggesting that the associated Cu 12e atom is r attling between two lone pair Sb. The number of locking and rattling Cu 12e a t different temperatures can also be found in Table 2 .5 . Locking 12e appeared when temperature went below MST temperature, and became more popular as temperature decrease. Notably, at 20 K, the distance between Cu 12e and lone pair Sb is ( 3.17 Å) larger th a n the other cases, which could compensate the shrinkage of unit cell and lead to negative thermal expansion. Table 2 . 5 : Number of rattling and locking Cu 12e, Cu 12e - Sb bond order and bond distance for Cu 12 S b 4 S 13 at different temperatures. Number of rattling Cu 12e Number of locking Cu 12e Locking Cu 12e - Sb bond order Cu 12e - lone pair Sb distance (Å) 20K 6 6 0.2169 3.17 50K 7 5 0.1784 3.06 87.5K 10 2 0.1386 3.07 300K 12 0 No locking 12e (lone pair) 3.12 40 (a) (b) (c) Figure 2 . 5 : Schematic of Cu - Sb interaction at 20 K and 300 K (a) and potential energy landscape based on Sb - Cu12e - Sb cluster a t (b) 20 K and (c) 300 K. 2. 3.3 Vibrational Density of States Calculated VDOS is derived from the Fourier transform of velocity autocorrelation function of atomic trajectories as following: , ( 22 ) where and are the mass and velocity of nth atom, the total number of atoms in the 41 group, and the Gaussian window function with a peak width of 1 meV 152 . Calculated partial vibrational density of states of each atomic group is presented in Figure 2 .6. It is obvious that the low energy vibr a tional mode (~3 meV at 300 K), which is contributed by Cu 12e atoms, is shifting to lower energy with decreasing tempera ture, i.e., phonon softening upon cooling 138 , 153 - 154 . The observation is consistent with experimental INS data 153 and validate our simulation results. Figure 2 . 6 : Partial vibration density of states (VDOS) of each atomic group for Cu 12 Sb 4 S 13 at different temperatures. 2. 3.4 Thermoelectric properties Density of states (DOS) calculation of Cu 12 Sb 4 S 13 at different temperatures, is pre sented in 42 Figure 2 .7 . We used the average structure of NVE run as input structure for each temperature. At 300 K, the Fermi level E F falls into top of valence bands with a sharp peak suggesting the degenerate d p - type semiconductor feature of Cu 12 Sb 4 S 13 with a bandgap of 1.30 eV, which is in good agreement with previous reports 35 , 40 , 155 - 156 . Overall feature do not vary a lot as temperature decreased, except that band gap energy slightly dropped from 1.30 eV to 1.25 eV. Figure 2 . 7 : DOS calculation for Cu 12 Sb 4 S 13 and partial contribution from each atomic group at different temperatures. Fermi energy is marked by dashed line. 43 Thermoelectric properties of Cu 12 Sb 4 S 13 are calculated by Boltzwann tran sport e quation implemented in VASP . In the system where an electric f ield and temperature gradient T exist, current density and the electronic heat current can be described as 157 - 158 T), ( 23 ) T, ( 24 ) where is electrical conductivity, is Seebeck coefficient and is the thermal conductivity. Since we only take the electronic contribution part account for , electrical thermal conductivity here is defined as heat current per unit of temperature gradient in open circuit condition and can be expressed as : , ( 25 ) Based on the semiclassic al transport theory using the Boltzmann transport equation presented by Pizzi et al 146 , the expre ssion for tensors , and can be described as a function of chemical potential and temperature : , ( 26 ) , ( 27 ) , ( 28 ) wh ere and are the Cartesian indice, the usual Fermi - Dirac distribution function and the transport distribution function (TDF). In Figure 2 .8 (a), the tempe r ature dependence of electrical resistivity has been showed. The calculated resistivity values from BoltzWann are roughly one orde r magnitude lower than the experimental data during the whole temperature range. Moreover, 44 study i s unable to catch, despite discontinuity can be observed at 87.5 K. This result is in good agreement with our previous DOS calcul ation. This is possibly due to the limited ability of using PBE as XC functional in this calculation. (a) (b) (c) - Figure 2 . 8 : Temperature dependencies of calculated (a) electrical resistivity, (b) Seebeck coefficient and (c) electrical thermal conduc t ivity for Cu 12 Sb 4 S 13 . Experimental data from Nasonova el al. 136 are shown here for comparison. 45 Seebeck coeffic ient S from BoltzWann calculation have been show n in in Figure 2.8(b). S maintained positive at the whole temperature range, suggesting that holes are major carriers in this structure. In addition, S increased as temperature increased above MST temperature , which is another evidence of degenerate p - type s emiconductor feature of Cu12Sb4S13. Comparing to experimental data, from MST temperature to 250 K, good agreement between simulation results and experimental values. Above 250 K, instead of keeping growing linearly, calculated data showed a decreased rat e of growing and approached a plateau at high temperature range. Below MST temperature, computed data did not have the local maxima at 67 K however the abrupt increasement at MST temperature can be observed. Thermal conductivity is composed of two parts, t he lattice part l and the electronic part e . In the experimental data, there was an upturn at MST temperature for e , which was contributed from the remarkable electrical conductivity of metal, and relatively smoother grow for l . Here we obtained e b y applying E quation ( 2 5 ) as shown in Figure 2 .8 (c). The upturn at MST temperature can be identified, however, the absolu te values were higher than the experimental data for the whole temperature range, due to the lower estimation of resistivity calculation . 2. 4 Conclusion In this work, we examined the low - temperature structural properties of pristine tetrahedrite Cu 12 Sb 4 S 13 by neutron diffraction measurement, along with DFT - based molecular dynamics simulation. No evidence of the reported cubic - tetragonal ph as e transition occurred accompanied with MST from our experimental and computational results. Negative thermal expansion below 50 K can be observed, which we ascribed to the interaction between locking Cu 12e and coordinated Sb atom. Density of states calcu la tions by DFT did not reveal MS T, but slightly bandgap reduction as temperature decreased can be observed. Moreover, we investigated the low 46 energy phonon vibrational mode which accounted for the anomalous phonon softening. We also calculated thermoelectr ic properties such as electrical resistivity, Seebeck coefficient and electronic thermal conductivity by BoltzWann code at low temperature and compared with experimental data. Further investigations are required to understand the origin of MST and the rela te d properties of Cu 12 Sb 4 S 13 . 47 CHAPTER 3 : Incoherent and coherent atomic dynamics of Cu 10.5 NiZn 0.5 Sb 4 S 13 3.1 Introduction Previous study indicated that the transition metal dopants such as Ni, Zn, Fe and Mn, can improve the thermoelectric performance o f tetrahedrites by adjusting el ectronic structure and reducing lattice thermal conductivity 35 - 36 . Previous experimental and computational studies found that the lo w thermal conductivity in tetrahedrites accompanied by anomalous atomic dynamics 45 , 152 . One of the anomalous behaviors is the unusually large anisotropic displacement (ADP) of trigonal planar Cu 12e atom s. Our previous study 152 indicated that the large ADP can be attributed to the Cu 12e - Sb lone pair. Weak bonding between Cu 12e and Sb leads to large vibration amplitude, results in the out - of - pla ne anharmonic rattling. The other type of anomalous behavior is called neutron scattering (INS) measurements. This behavior can be described as low - energy vibra tion mode energies decreased with the decreasing temperature. Although it has been argued that this phenomenon is also due to the dynamics of Cu 12e atoms, the details remain unknown. This work will apply DFT - based MD simulation, along with inelastic neutro n scattering ( INS) experiment, to study the atomic mechanism of Cu 10. 5 NiZn 0.5 Sb 4 S 13 tetrahedrites 154 . 3.2 Compu tati onal and experimental details Ni - doped tetrahedrite with the composition Cu 10. 5 NiZn 0.5 Sb 4 S 13 has been chosen as a 34 . DF T - based MD simulation was performed with t he Vienne Ab initio Simulation Package (VASP ) 143 - 144 , 159 - 160 employing the Projector Augmented - Wave (PAW) method 161 - 162 . GGA of PBE parametrization 100 was used for exchange - correlation functional. Simulation cell was 1 x 1 x 1 supercell with 58 atoms in total, while Ni and Zn are randomly distributed at Cu 12d sties. 48 Valence configurations are 3p 6 4s 1 3d 10 , 4s 1 3d 9 , 4s 2 3d 10 , 5s 2 5p 3 a nd 3s 2 3p 4 fo r Cu, Ni, Zn, Sb and S respectively. First a constant number of particles, pressure and temperature (NPT) ensemble was implemented for 3 ps, with a Langevin thermostat and Parrinello - Rahman barostat, to obtain lattice parameters at various temp erature. Nex t a constant number, volume, and energy (NVE) ensemble was performed for 5 ps for equilibration and followed by 20 ps production run for sampling atomic trajectories. Plane wave energy cutoff was 450 eV, with a gamma point K - mesh, and timestep was 1 fs. IN S of a powdered sample at 300 K were implemented by the wide Angular - Range Chopper Spectrometer (ARCS) at the Spallation Neutron Source of Oak Ridge National Laboratory, with incident neutron energy of 20 meV. 3.3 Results and discussion 3 .3.1 I ncoherent at omic dynamics of Cu 10.5 NiZn 0.5 Sb 4 S 13 The calculated vibrational density of states (VDOS) is calculated from MD trajectory using the method as shown in Section 2.3.3, to compare the one extracted from INS measurement, as shown in Figure 3 .1(a). VDOS from IN S measurement is obtained by making the one - phonon incoherent scattering approximation 163 , with Debye - Waller factors from MD trajectories. The comparison indicates relatively reasonable agreement between computation and ex periment, while the high energy portions difference is due to the one - phonon incoherent ap proximation. Figure 3 .1 ( b) shows the partial VDOS spectra for each atom group. It is worth noting that a low energy peak exhibited at around 4 meV, which is predomina tely contributed by Cu 12e atoms. This low energy vibration mode often suggest s weak bondin g and large - amplitude, and possibly suppresses the acoustic phonon branch, and leads to low thermal conductivity. 49 Figure 3 . 1 : (a)Comparison of VDOS from MD simulation and INS experiment. (b) P artial VDOS f or each atom group and total VDOS from MD simulation at 300K. 50 Figure 3 . 2 : Partial VDOS of Cu 12 e atoms at different temperatures. Shape was fitted to sum of two Lorentzian functions. Previous e x perimental IN S data indicated that the low energy modes shifted to lower energies with lower temperature, i.e., phonon softening upon cooling, for Cu 10 Zn 2 Sb 4 S 13 and Cu 10 Zn 2 Sb 2 Te 2 S 13 138 , 153 . To understand the mechanisms, we calculate p artial VDOS f or Cu 12e ato ms at different temperatures, as shown in Figure 3 .2. Two peaks (at ~3.25 eV and ~4.8 eV ) can be identified for 50 K, 100 K and 200 K, which can be fitted to the sum of two Lorentzian peaks. The first peak located at 3.25 eV for 50 K and 100 K, then shifted t o higher energy as temperature increases, which is consistent with experimental INS d ata. Our previous work 152 suggested that the Sb lone pair electrons modulate the motion of Cu 12e atom and lea d to its anharmonic rattling. To visualize to interaction between Cu 12e and Sb, we co nvert the distribution of distances to effective one - pair - potential (Figure 3 .3(a) and 3 .3(c)) and correlated to VDOS of Cu 12e (Figure 3 .3(b) and 3 .3(d)). At 50 K, a doubl e - well shape (highlighted red in Figure 3 .3 (a)) can be identified 51 for the Cu and Sb interaction energy. This double - well shape is consistent with previous lattice dynamic calculation 164 , w hich indicates that the associated Cu 12e atom is locked to one of the bonded Sb with average distance of 2.9 Å (Figure 3 .3(e)). The Density Derived Electr ostatic and Chemi cal (DDEC) scheme 151 calculation show that Cu 12e - Sb bond order is 0.25, more like a bond ing pair instead of lone pair. The certain Cu 12e thus contributed to high vibrational energy of 5.5 meV as Figure 3 artic shape (highl ighted blue in Figure 3 .3 (b)) is also observed, which suggests that associated Cu 12e is rattling between two lone pair Sb with average bond order of 0.1, and is responsible for the low vibrational energy of 3.25 meV. Then correlation at 300 K has been st udied as Figure 3 .3(c) and Figure 3 .3(d). All Cu 12e atoms at 300 K exhibit as rattler s and vibrational amplitude increased by additional thermal energy, thus results in higher vibrational energy, i.e., softening upon cooling. Figure 3 . 3 : Potential energy landscape based on Sb - Cu 12e - Sb clus ter at 50 K (a) and 300 K (b), and correlated VDOS spectra of Cu 12e at 50 K (c) and 300 K (d). Cu 12e - Sb interaction are shown in (e). 52 3 . 3 .2 Coherent atom ic dynamics of Cu 1 0.5 NiZn 0.5 Sb 4 S 13 Figure 3 . 4 : Calculated dynamic structure factors based on the momentum correlation function (a) longitudinal and (b) transverse direction. Coherent atomic dynamics are a lso examined by co mputing the longitudinal and transverse dynamics structure factor based on the momentum correlation function, as shown in Figure 3 .4. The minimum Q value in this work is 0.6 Å - 1 due to the limited cell size. A clear optic band can be iden tified at ~ 4 meV. Because of the acoustic - optic interaction, i.e. 165 - 166 , 53 dispersive bands at ~10 - 20 meV/0.6 - 1.3 Å - 1 for longitudinal, ~5 - 10 meV/0.6 - 1.8 Å - 1 for transverse direction, are rem iniscent of the opt is not r ecognizable, which suggests that it is confined in a small dynamic region which is less than 4 meV and 0.6 Å - 1 on mean free path i s close to the interatomic spacing and lead to low thermal conductivity. To validate the significance of the hypothesis s emi - quantitatively, we obtaine a Debye temperature of 281 K and Debye heat capacity of 2.81 k B /atom at 300K. Althoug h our Debye tempera ture is much higher than the value of 500 K from heat capacity studies 167 , our Debye heat capacity is comparable to the value of 2.75 k B /atom. I n addition, the mean sound speed have been computed with the value of 2522 m/s, which is close to the sound speed of 2668 m/s in hydrodynamic region from bulk modulus 168 . Moreover, by making the approximation tha t the mean free path equals to the average interatomic spacing, we can obtain the thermal conductivity of 0.45 W m - 1 K - 1 , which is simi lar to the experimental value of 0.43 W m - 1 K - 1 34 . The good agreements between our study and literatures provide the validation of previous hypothesis. 3.4 Summary In t his work we investigate the incoherent and cohere nt atomic dynamics of Cu 10.5 NiZn 0.5 Sb 4 S 13 tetrahedrite using first - principle molecular dynamics simulation. Fo r the incoherent dynamics, the computed vibrational density of states is in good agreement with the experimental result from inelastic neutron dif fraction data at 300 K. Detailed analysis of Cu motion inside the Sb[CuS 3 ]Cu atomic cage indicates the reduct ion of anharmonic rattling of Cu between two Sb atoms upon cooling, which leads to anomalous phonon softening. For the coher ent dynamics, the dyna mic structure factors in the longitudinal and transverse direction suggest that the acoustic modes are confine d in a small region of the scattering space, which accounts for the 54 small thermal conductivity measured. This study could fu rther our understandin g of the atomic dynamics of tetrahedrite thermoelectrics and can more generally provide insight on how materia ls can be designed at the atomic level to possess intrinsically low lattice thermal conductivity, an important parameter for not only thermoelectr ics but also other applications such as thermal barrier coatings. Th is work was published on APL Materials 4. 10 (2016) :104811. https://doi.org/10.1063/1.4959961 55 CHAPTER 4 : Mo bile Cu moveme nt in Cu - rich tetrahed rite Cu 14 Sb 4 S 13 4.1 Introduction The Cu - Sb - Cu phase diagram indicated that tetrahedrite composition varies from Cu 12 Sb 4 S 13 to Cu 14 Sb 4 S 13 depending on temperature 169 . Tetrahedrites with more than 12 Cu, i.e., Cu - rich tetrahedrite, had better thermoelectric performance than the pristine cell according to previous study by Vaqueiro et al 170 . The thermoelectric figure of merit for Cu 14 Sb 4 S 13 reached ~ 0.6 at 573 K, wh ich was similar to reported substituted tetrahedrites Cu 11 Zn 1 Sb 4 S 13 164 and Cu 12 Sb 3 Te S 13 41 . Vaqueiro also s uggested that, for Cu - rich tetrahedrites, Cu becomes mobile above 393 K and are likely to be su perionic conductors, while the excess Cu atoms were likely to be at Cu 24 g sites. However , the insight of how Cu atoms distribute and move remains unknown. Inspired by the previous work, it is interesting to investigate the details of mobile Cu movemen t, wi th the help of DFT - based MD simulation using a Cu - rich tetrahedrite Cu 14 Sb 4 S 13 as a model . To help us understan ding the local environment of different Cu sites and possible Cu diffusion pathways, we manually put extra Cu atoms into a stoichiometric Cu 1 2 Sb 4 S 13 cell (I - 43m) and make Cu 24g sites fully occupied with the composition of Cu 1 4 Sb 4 S 13 , using a three - dimensional visualization software VESTA 37 . The initial crystal structure and representative unit around Cu 24g are shown in Figure 4 .1 and Figure 4 .2. Three distinct Cu sites can be identified, as Cu12d, Cu 12e and Cu 24g site. The local environment of Cu 12d and C u 12e have been discussed in Section 1.2.3 , while the additional Cu 24g are bonded to three S atoms and encompassed with two C u 12d and t wo Cu 12e atoms. From the crystal structure of Cu 1 4 Sb 4 S 13 , it is reasonable to hypothesize that mobile Cu atom may move ins i de the structure by jumping between 12d , 12e and 24g sites. 56 Figure 4 . 1 : Initial crystal structure of Cu - rich tetrahedrite Cu 14 Sb 4 S 13 with three distinct Cu sites. Cu24g atoms are displayed as black. F i gure 4 . 2 : Cu24g local environment for better visualization. 4.2 Computational details Cu - rich tetrahedrite Cu 14 Sb 4 S 13 has been chosen as the model structure in the study, and excess Cu atoms are distributed i n Cu 24g sites based on previous research 155 , 170 - 171 . Simulation setups are the same as Section 3 .2, exc ept that PBE - D3 172 is used as exchange - correlation functional with 2x2x2 K - po ints mesh. NVE pro duction simulation time is increased to 40 ps in 57 this case. Simulation temperature is chosen to be 700 K, since Cu - rich tetrahedrite will collapse into a single phase above ~ 500 K based on the neutron diffraction data 170 . 4.3 Results and discussion 4.3.1 Lattice parameters Lattice parameter obtained f rom N PT simulation is 10.3972 Å at 700 K, comparing to 10.4924 Å at 673 K from literature. Regarding that PBE - D3 method usually underestimates lattice parameters for about 4%, the difference is considered to be acceptable 173 . Lattice fluctuation during NPT can be observed as Figure 4 .3 . The relatively large fl uctuat ion suggested the structural instability, which we ascribe to the mobile Cu atoms movement inside the cell. Figure 4 . 3 : Lattice parameters (a, b, c) as function of time in NPT simulation. 58 4.3.2 Incoh erent density correlation and Cu self - diffusivity Here we use incoherent density correlation ( ) to investigate the diffusivity of Cu - rich tetrahedrite. is calculated from MD trajectory for each Q value 107 , with correlation time of 35 ps as following: . ( 29 ) where is the atomic specie, the number of atoms, the wave vector and the position of n th atom at time t. of diffusing species will exhibit a c lear de cay to 0 (ideally), while the vibrating ones will show rapid decay and then non - zero plateau. Simulation cell have 28 Cu atoms in total. After examining the shape of calculated for each Cu atom, these Cu atoms can be divided in to two group: vibrating Cu and di ffusing Cu. Then of all vibrating Cu are summed and averaged, as well as the diffusing group, as shown in Figure 4 .4, where only three smallest values are inlcuded. Ignoring the first 5 ps, where at oms are undergoing ballistic collision and other short - time dynamics, the decay of diffusing Cu group plot can be described as a stretched exponential function 174 : , ( 30 ) where is scaling factor, 1/ is the relaxation time and is stretching parameter. vs. 2 plot is shown in Figure 4 .5 , a linear dependence of with ( = ) is characteristic of continuous random walk translational diffusion, i.e., Fickian diffustion, where the slope D i s self - diffusion constant 107 . In addition to applying Fickian diffusion model to describe the continuous random walk translational diffusion at low Q range, the broadening of with Q 2 reache d a plateau asymptotically at higher Q , i n dicating a typical jump - diffusion process, and can be described as 59 the Singwi - Sjölander (SS) jump model 175 - 176 in which diffusion occurs via jump diffusi on where ( 31 ) where is the residence time and the mean square jump dist ance . At low value SS model reduces to Fickian model with and jump length . Cu diffusivity obtained from Fickian model and SS model yield 4.8961 x 10 - 6 cm 2 /s and 1.0 x 10 - 5 cm 2 /s, r espectively. Residence time and mean jump l ength can be extracted by SS model fit, with values of 0.9885 ps and 0.9947 Å, resepectively. (a) (b) Figure 4 . 4 : Incoherent density correlation for (a) diffusing Cu atoms and (b) vibrating Cu atoms f or three small Q values. 60 Figure 4 . 5 : HWHM ( ) as a function of . Red line is fitted to Fickian model and b lue line is fitted to Singwi - Sjolander model. 4.3.3 Cu nuclear density map Although the quantitatively incoherent density correlation analysis offered useful information including Cu sel f - diffusivity, physical det ails such as the distribution and movement of Cu atoms remain unknown. Nuclear density map enables the visualization of atom distribution by dividing the cell into cubic pixels of roughly 0.1 Å in length, then calculate the proba bility density of each type of atom in each pixel 177 . Figure 4 .6 was generated by VESTA, showing the density maps of representative units of diffusing pathways. Isosurface level is determined to be 0.12 Å - 3 f o r better visualization. From the nuclear density map, various diffusion paths for Cu atoms can be identified. For example, a diffusion path between 12d site and 24g site can be observed as Figure 4 .6(a), the intermediate site could be a 48g site. We also n oticed that Cu atom can jump from a 12d cage to another 12d cage, without taking a 12e or 24g site as intermediate 61 transition (Figure 4 .6(b)). In addition, Cu atom can diffuse between 12e site and 24g site back and forth (Figure 4 .6(c)). No nuclear densi t y was identified between Cu 12d and Cu 12e sites, which means a direct jump between these two sites is not favorable for Cu atoms in this case. It is noticeable that one of the Cu 12d sites is never accessed by Cu atom during the whole simulation time (40 ps i n total). We think that might due to the finite transport ability of Cu at 700 K, or limited statistics by using 1x1x1 simulation cell and 40 ps simulation time. Moreover, identified diffusion paths with the participation of atoms in Cu 12d and Cu 24g site s are much more than the paths including Cu 12e atoms, suggesting that Cu atoms are stable at Cu 12e sites, which is interesting considering the weak bonding of Cu 12e Sb inside Sb[CuS 3 ]Sb cage, and need to be further investigated. (a) (b) (c) Figure 4 . 6 : Nuclear density may of Cu atom showing the diffusion pathway: (a) 12d 24g path; (b) 12d 12d path; (c) 12e 24g path. Brown atom is Sb and yellow is S. Isosurface level of 0.12 Å - 3 . 62 4.3.4 Energy barriers Figure 4 . 7 : Diffusion paths schematic of (a) 24g 12d and (b) 24g 12e extracted from the actual MD trajectory and m igration energy barriers of (c) 24g 12d and (d) 24g 12e from NEB calculations. In order to estimate the migration energy barriers of different migration paths, here we NEB ) 178 - 180 methods as programmed in the VASP package, using 5 images between each metastable configuration. Two diffusion pat hs, 24g 12d and 24g 12e 63 are selected as examples for comparison. The initial sta te s and final states are chosen among our MD trajectory, instead of manually creating vacancies, which leads to better estimation of the real case. Zero of energy is chosen as the initial state energy. Figure 4 .7 (a) and (b), obtained by CrystalMaker 181 , re veal the schematic of those two diffusion paths, where only part of the crystal cell are shown for better visulization. Figure 4 .7 (c) displays the migration energy diagram of 24g 12d , with migration energy barrier E m of 1.32 eV, where two intermediat e l ocal minimums (corresponding to Cu 48h sties) can be observed. Whereas, E m of 24g 12 e is estimated to be 1.56 eV as displayed in Figure 4 .7 (d), which is higher than the 24g 12d case, indicating the 24g 12 d diffusion pathway is energetically prefer abl e than 24g 12 e . This result is in g ood agreement with the above density map analysis that less diffusion paths involving 12 e atoms can be observed. 4.4 Summary In this chapter we discussed the Cu ion transport properties of copper - rich tetrahedrite Cu 1 4 S b 4 S 13 with the help of DFT - based MD. Large fluctuation of lattice parameter from NPT indicated the structural instability during the simulation. By investigating the shape of incoherent density correlation at 700 K, we can conclude th at some Cu atoms remain vibrating while the other become mobile by diffusing among different Cu sites. Diffusivity, residence time and mean square jump length of mobile Cu atoms can be extract by fitting simulation data to Fickian model and SS model. Th e n uclear den sity map is also examined to visualize the motion of Cu atoms at different sites. Cu 12d and Cu 24g tend to be more mobile than Cu in 12e sites within the simulation. Future study with larger cell and longer simulation time will be needed to con fir m this. Mo reover, several diffusion paths can be identified, including direct jump between two 12d E m for 24g 12d and 24g 12e diffusion paths, with t he values of 1.32 eV and 1.56 eV, 64 respectively. Further study of nonstoichiometric tetrahedrite as potential ionic conductor will be interesting, including optimization of Cu - rich composition and experimental validation of the trans port properties such as Cu self - diffusivity, ionic conductivity and activation energy. 65 CHAPTER 5 : Density - f unctional t heory b ased m olecular d ynamics s imulation of t etrahedrite t hermoelectrics: effect of cell size and basis sets 5.1 Introduction With the a dvancement of hardware and sof twa re infrastructure, computational modeling has become an indispensable tool in materials research, often in parallel with experimental effort. Among different modeling methods, density - functional theory (DFT) in the Kohn - Sh am approach 182 has been widely used in the investigation of chemicals, hard materials, and soft matter. Although DFT calculations are considered as first - principles methods they are not optio n fr ee, as the size of simulation cell, basis sets for the wavefunction, excha nge - correlation (XC) functional, etc. must be selected first. While there are many studies on the effect of XC functionals, e.g. the recent study by Tran et al. 183 , the ex amina tion of cell size is rare due to the computational overhead associated wit h the increase of atoms/electrons. One of the size studies by Spiekermann et al. 184 suggested different trends, e.g. water release for a 192 - atom cell and water uptake for a 96 - atom cell of supercritical H 2 O - SiO 2 fluids, pointing to a finite size effect. In terms of basis sets, plane waves (PW) and atomic orbitals (AO) are two popular choices. Recently, Miceli et al. 185 compared the structural, dynamic, and electronic properties of liquid water using both AO and PW basis sets and found good agreement on results from two sets. Ulian et al. 186 also compared PW and AO basis sets using Mg 3 Si 4 O 10 (OH) 2 layer silicate and found both basis sets a dequately describe the geometry, energy , and infrared spectra. Inspired by these studies, in this work we report our investigation of effect of cell size and basis set on re sults fr om DFT - based first - principles molecular dynamics simulation of a thermoele ctric material: Cu 10 Zn 2 Sb 4 S 13 tetrahedrite. Tetrahedrite materials with a general composition of Cu 12 - x M x Sb 4 S 13 where M is a metal dopant have emerged recently as promising 66 thermoe lectric candidates due to their elemental abundance, environmental friendliness, favorable electronic properties, and most importantly low lattice thermal conductivity (< 1 W m - 1 K - 1 for a wide temperature range) 164 , 187 . While this group of mat erials has been the subject of several DFT calculations 34 , 42 , 164 , 188 - 192 cell size and basis set remains elusive. We selected two cell sizes (1x1x1 and 2x2x2) and two basis sets (PW and AO) and compared simulation results from three simulations (1x1x1 PW, 1x1x1 AO, and 2x2x2 AO) and with experiments. 5.2 Computational and exp erimental deta ils Vienna Ab initio Simulation Package ( VASP ) 193 - 196 code employing the Projector Augmented - Wave (PAW) method 197 - 198 and Quickstep code 199 implemented in cp2k 200 are employed as DFT pac kages using PW and AO basis sets, respectively. In VASP , valence electron configurations for Cu, Zn, Sb, S atoms are 4s 1 3d 10 , 4s 2 3d 10 , 5s 2 5p 3 , and 3s 2 3p 4 , respectively. The plane wave energy cutoff was 450 eV. In Quickstep, the mixed Gaussian and plane w ave (GPW) method, i.e. Gaussian basis sets for orbitals with auxiliary plane waves for electr on density, was used. Valence electron configurations were the same as in VASP but with Godecker Tetter Hutter (GTH) norm - conserving pseudo potentials 201 - 202 . The plane wave cutoff was 1200 Ry and the atomic orbital (Gaussian) basis sets were molecular optimized double zeta - valence basis sets with a polarization fun ction (DZVP) 203 . In bot h cases, the generaliz ed gradient approximation (GGA) with Perdew - Burke - Ernzerhof (PBE) parametrization 204 , 205 was used for the XC functional and a single Gamma point was used for the k - mesh. Two sizes of simulati on cells, 1x1x1 (58 atoms) and 2x2x2 (464 atoms) were used based on the crystal structure of tetrahedrite (I - 43m) with Zn randomly distributed at Cu 12d sites. We chose 1x1x1 cell to compare AO and PW basis sets and 67 AO basis set to co mpare 1x1x1 and 2x2x2 cells. Due to the exorbitant computational cost, we did not use the 2x2x2 cell and PW combination. For MD simulation, first a constant number of particles, pressure and temperature (NPT) ensemble at 300 K and zero pressure was implem ented for 3 ps for 1x1 x1 cells in order to obtain the lattice parameters. The Langevin thermostat (friction coefficient of 50 and 1 ps - 1 for atoms and lattice) and Parrinello and Rahman barostat (lattice mass of 10 a.m.u.) was used in VASP , while a Nose t hermostat (time consta nt of 0.1 ps) and Martyna - Tuckerman - Tobias - Klein (MTTK) barostat 206 (time constant of 0.5 ps) was used in Quickstep. Next a constant number, volume, energy (NVE) ensemble was performed for 5 ps for equilibration and 20 ps for production at 300 K, while in NVE there is no temperature cor rection was made. The time step was 1 fs for both NPT and NVE runs. Synthesis of Cu 10 Zn 2 Sb 4 S 13 powder was performed by mechanical alloying as described in the literature 51 . Cu K - edge X - ray absorption exp eriments were performed at the 4 - 3 beamline of Stanford Synchrotron Radiation Laboratory. We used the software SIXPack 207 to average raw data files out of three absorption scans and employed the software Athena 208 to extract the EXAFS signal. 5.3 Results and discussion 5.3.1 Average structure Average structure of Cu 10 Zn 2 Sb 4 S 13 after 1x1x1 PW NVE simulation is shown in Figure 5 . 1. Two distinct C u sites can be identified, as Cu 12d and Cu 12e . Cu 12d atom coordinated by four S atoms and form a tetrahedron, while Cu 12e atom surrounded by three S atoms and form a triangular planar coo rdination. All the Zn dopants located in Cu 12d sites. 68 Figure 5 . 1 : Average structure after 1x1x1 PW NVE simulation showing the crystal structure of Cu 10 Zn 2 Sb 4 S 13 . 5.3.2 Lattice parameters (a) (b) Figure 5 . 2 : Lattice para meters ( a , b , c ) of (a) 1x1x1 PW and (b) 1x1x1 AO cell as a function of time in NPT simulation. 69 Lattice parameters ( a , b , c ) from the NPT ensemble of 1x1x1 PW and 1x1x1 AO simulations at 300 K are shown in Figure 3 . 1 . It can be seen both simulations suggest a cubic structure. The average cubic lattice parameters are 10.4692 0.135 Å and 10.5113 0.073 Å, respectively. They are close to each other within the standard deviation. They ar e both higher than the lit erature value of 10.3833 Å 209 , as the PBE X C functional tends to overestimate lattice parameters for various materials 183 . We took the lattice parameter of 2x2x2 AO cell as double of 1x1x1 AO cell. 5.3.3 Vibrational density of states (a) (b) (c) (d) (e) (f) Figure 5 . 3 : Partial and total vibrational density of states (VDOS) of each atom group (a - e: Cu 12e , Cu 12d , Zn 12d , Sb, and S; f: total) for three simulations (1x1x1 PW, 1x1x1 AO, and 2x2x2 AO) . We inv estigated the effect of cell size and basis set on vibrational density of states (VDOS), 70 which can be extracted from the atomic trajectory as discussed in Section 2.3.3. The partial VDOS (a - e) of each atom group, Cu 12e , Cu 12d , Zn 12d , Sb, and S and total VD OS (f) are shown in Figure 5 . 3 . The 12e and 12d sites are trigonal planar and tetrahedral positions, respectively. In Figure 5 . 3 ( a ) , the vibration of Cu 12e shows three pe aks (~4.2, 5.9, 7.5 meV) at the low - energy portion with the first peak being strongest for 1x1x1 PW simulation, while 1x1x1 AO simulation has three peaks (~4.8, 6.1, 7.8 meV) with the second and third peaks almost equal in intensity. If we compare 1x1x1 A O and 2x2x2 AO simulations, we can see the peak at 7.8 meV from the former decays into a shoulder in the latter. The similar observation was made for other atom groups (b - e) and the total group (f). We think that different basis sets cause slightly diffe rent electronic structures a nd atomic forces, leading to slight variation in vibrational properties. On the other hand, larger cells allow more atoms and vibrational modes to be sampled, which effect ively broaden the peaks. Under the harmonic approximatio n, the vibrational energy is connected to the total VDOS as the following: , ( 32 ) from which heat capacity can be computed. The constant - volume heat capacity per atom for three simulati ons is shown in Figure 5 . 4 . Heat capacities of t hree simulations are close to each other at all temperatures and approach 3 k B at high temperature according to the Dulong - Petit law. In the meantime, experimental da ta for the same compound measured by Lara - Curzio et al . 210 and for a similar composition Cu 10.5 Zn 1.5 Sb 4 S 13 measured by Lu et al. 164 are also included in Figure 5 .4 . Experimental values agree well with computed ones. It is worth noting that we assume the total VDOS ( obtained from 300 K, i. e. Figure 5 .3( f)) stay the same for the entire computed range, 0 71 900 K. Figure 5 . 4 : Heat capacity per atom for three simulations (1x1x1 PW, 1x1x1 AO, and 2x2x2 AO) . Experimenta l data from Lu et al. 7 and Lara - Curzio et al. 33 are also shown for comparison. 5.3.4 Phonon dispersion In addition to VDOS, we also examined the effect of cell size and basis set on phonon dispersion by comparing the Q - resolved coherent velocity correlatio n. For the coherent v elocity , ( 33 ) where is the total number of atoms, the longitudinal (L) and transverse (T) correlation can be calculated as the following: , ( 34 ) . ( 35 ) The allowed vectors are based on the simulation cell as where are the 72 reciprocal l attice vectors and are integers. While the larger cell size of 2x2x2 AO simulation enables smaller Q value to be studied, we picked two smallest values that are common to all three simulations (0.60 and 0.84 Å - 1 ) and plotted the longitudinal and transv erse response in Figure 5 . 5 . Similar to what we observed for VDOS in Figure 5 . 3 , results from three simulations are comparable but not identical. Examples of notable devi ations include the large intensity difference of transverse peak at ~2 THz in Figure 5 . 5 (b ) and of longitudinal peak at ~40 meV in Figure 5 . 5 (c). (a) longitudinal (b) tra nsverse (c) longitudinal (d) transverse Figure 5 . 5 : Longitudinal phonon dispersion for (a) Q=0.60 A - 1 , (c) Q=0.84 A - 1 and transverse phonon dispersion for (b) Q=0.60 A - 1 , (d) Q=0.84 A - 1 . 5.3.5 EX AFS spectra Extended X - ray Absorption Fine Structure (EXAFS) experiments are often applied to study the local structure around an atomic species. In this work we will also compare Cu edge EXAFS 73 spectra from computation and experiment. Computationally, th e EXAFS spec trum was obtained from configuration averaging of atomic trajectory. Briefly, we selected 1 frame for every 10 fs of the trajectory to have a total of 2,000 frames. In each frame, we defined an atomic cluster around each Cu atom (20 in total) with a 5 Å r adius and employed the software FEFF9 211 to calculate the EXAFS spectrum. We averaged over all selected frames and all Cu atoms in each frame and plotted the result in Figure 5 . 6 . Similar to the heat capacity p lot in Figure 5 . 4 , results from three simulations are close to each other. The peak position difference between computation and experiment in Figure 5 . 6 is likely caused b y the difference in lattice parameter, i.e. Section 5 . 3. 2 , as these peak positions are related to specific bond lengths. Figure 5 . 6 : k 2 - weighted Cu EXAFS spectrum obtained by experiment (black circle) and s imulations (solid lines). 5.3.6 2x2x2 AO Simulation The low Q values accessible only to the 2x2x2 cell allow us to examine the region closer 74 to the h ydrodynamic limit of atomic dynamics, where useful materials properties could be extracted 212 . First, we could study sound speeds by examining the longitudinal and transverse momentum correlation. They have the same expression as in Eq. (27) but with the expression as : . ( 36 ) Longitudinal and transverse response of coherent momentum correlation with three smallest Q values corresponding to the 2x2x2 cell are shown in Figure 5 . 7 . In the longitudinal case, the strongest peak of each curve is moving to higher energy with increasing Q . We hypothesize that this peak corresponds to the longitudinal sound mode and plot i ts energy as a function of Q in Figure 5 . 7 ( c ) . A linear fit through the origin yields a longitudinal sound speed of 3945 m/s. For the transv erse response, two peaks can be observed and we ascribe them to transverse sound modes in Figure 5 . 7 ( b ) , with the Q dependence of peak position shown in Figure 5 . 7 ( c ) . It is possible that the two smallest Q are in the linear dispersion region but we only used slopes between the origin and smallest Q to obtain sound speeds of 1805 and 2186 m/s. We took the average of these two, 1995 m/s , as the transverse sound speed . Sound speeds in cubic crystals are related to elastic properties su ch as shear ( G ) and bulk ( B ) modulus as the following: ˆx ˆx , ( 37 ) where is the mass density, and the elastic constants . The density from our 2x2x2 AO simulation is 4.77 g/cm 3 , while literature values range from 4.67 to 5.05 g/cm 3 due to sample variation 168 . Our calculated shear and bulk modulus are 1 9 .0 GPa and 48.9 GPa , respectively, which are similar to experimental values of 19.1 23.0 GPa and 46.3 56.3 GPa 168 . 75 (a) (b) (c) Figure 5 . 7 : Longitudinal (a) and transverse (b) response of the momentum correlation. (c) Dispersion of longitudinal (L) and transverse (T) sound modes. Peak energy for each Q value (red dots: longitudinal; blue dots: transverse) and their linear regression have been shown. Second, this larger cell also enables t he id entification of acoustic modes in phonon dispersion. Similar to Figure 5 . 5 , longitudinal and transverse response of coherent velocity correlation with three smallest Q values corresponding to the 2x2x2 cell ar e shown in Figure 5 . 8 . It is immediately recognizable that coherent momentum ( Figure 5 . 7 ) and coherent velocity ( Figure 5 . 8 ) re sp onses are similar at these Q values. This is not surprising as it is expected that mass variation from different species in a multicomponent material will only have small influence at t he hydrodynamic limit, i.e. the material behaving as a single compon en t material with effective mass. The closeness of elastic properties between experiment and computation leads us to assign the strongest peak in Figure 5 . 8 ( a ) as the longitudinal acoustic (LA) mode. In addition, p ositions of two shoulder peaks ahead of the LA peak, ~ 4 and 6 meV , do not vary significantly with Q so we described them as quasi - localized (QL) modes with optic character. Such pre - LA peaks show up as Boson peaks in VDOS (e.g. Figure 5 . 3 ( f ) ) and a comparison to Figure 5 . 8 ( b ) indicates that they are coupled to transverse acoustic modes, which is consistent with t he investigation by Shintani and Tanaka 213 in revealing a universal link between boson peak and transverse phonons in glass. 76 In the context of thermal properties, we belie ve these quasi - localized modes scatter heat - carrying acoustic modes by reducing their lifetime (peak broadening in Figure 5 . 8 ) instead of reducing their speeds, leading to a small thermal diffusivity and lattice th erm al conductivity. (a) (b) Figure 5 . 8 : Longitudinal (a) and transverse (b) response of the coherent velocity correlation. Longitudinal acoustic (LA) and quasi - localized (QL) modes are marked in a. T he dotted lines have the same energy betw een a and b. If we compare Figure 5 . 5 5.5 , Figure 5 . 7 , and Figure 5 . 8 , we can appreciate t he advantage of employing a larger simulation cell, i.e. allowing smaller Q value to be accessed, where it becomes easier to identify the acoustic/sound modes. For exampl e, the linear dispersion in Figure 5 . 7 coul d n ot to be observed in 1x1x1 cells. 5.4 Summary In this work we investigated the effect of simulation cell size and basis sets on the DFT - based molecular dynamics simulation results, using tetrahedrite Cu 10 Zn 2 Sb 4 S 13 thermoelectric as a model material. Th ree simulations, 1x1x1 cell with plane wave (PW) basis, 1x1x1 cell with 77 atomic orbital (AO) basis, and 2x2x2 cell with AO basis, were performed at 300 K. While variations can be observed for different cell sizes and basis sets, results of various structur al and dynamic properties, such as lattice parameters, partial and total vibrational density of states, phonon dispersion, heat capacity, and EXAFS spectra are close to each other for three simulations. The larger 2 x2x2 cell allows smaller Q to be accesse d ( closer to the hydrodynamic limit), where longitudinal and transverse sound/acoustic modes with linear dispersion can be identified. The bulk and shear modulus from extracted sound speeds agree with the experiment . In addition, two low - energy quasi - loc ali zed vibrational modes were detected. We believe that they scatter heat - carrying acoustic modes by reducing their lifetime instead of reducing their speeds, leading to a small thermal diffusivity and lattice therm al conductivity. This work was publishe d o n Computational Materials Science 144 (2018): 315 - 321. https://doi.org/10.1016/j.commatsci.2017.12.047 78 CHAPTER 6 : Structure and i onic c onduction s tudy on Li 3 PO 4 and LiPON with density functi ona l tight binding (DFTB) method 6 .1 Introduction Despite that no conclusive explanation has been made, N in the LiPON is commonly ascribed as the key to the high ionic conductivity and outstanding el ectrochemical stability. Researches have shown that the co nductivity of LiPON films increased as the N content increased, from 7 x 10 - 8 S/cm without any N content to 3.3 x 10 - 6 S/cm with 6% of N content 91 , 1 40 . However, how the N atoms coordinated with P atoms remains uncertain 214 - 217 . By analyzing N1s X - ray photoelectron spectroscopy (XPS) data on LiPON fi lms, B ates et al. indicated that N bridged among P atoms with double (N d ) or triple (N t ) coordination, and these crosslinked structure might cause the increase in Li + mobility by providing interconnected pathways 78 , 140 . However, another argument was made by Wang et al. that most N would form PO 3 N with apical N (N a ) and the rest would be N d 218 . Recent computational study supported the argument that only N a and N d configurations existed, while the rati o between this two was dependent on the Li content, with optimal composition Li 2.94 PO 3.50 N 0.31 in terms of ionic conductivity 219 - 220 . Investigating the both structures of the precursor Li 3 PO 4 and LiPON is the key to improve the battery performance and help illuminating the direction of future battery design. Density functional theory 221 (DFT) based molecular dynamics (MD) simulations are commonly used in accessing the structural, electronic and transport properties of materials. A lot of researche s on LiPON electrolyte have been done by DFT to simulate the crystalline (c - ) /amorphous (a - ) structures 219 - 220 , 222 , energy calculations 22 3 and the in terfacial reaction between LiPON and electrodes 224 - 225 . To study the transport properties, such as diffusion and ionic conduction, a large cell and a longtime scale simu lation (more samplings) are preferable to obtain 79 reliable results. Previous study in Chapter 5 also showed the adv antages of having a larger simulation cell size. Although DFT is known as accurate, it is limited by system size and simulation time. Thus, a suitable method which can balance the accuracy and efficiency is necessary. For example, here we employ a self - cha rge consistent density functional tight - binding (SCC - DFTB) 113 approach to study the structural and tran sport properties of Li 3 PO 4 and LiPON . 6 . 2 Computational details Details of DFTB theories and equations derivation c an be found in the introduction chapter, as well as at various publications 113 , 117 , 226 and will not be included here . All SCC - DFTB simulations were carried out by DFTB+ code proposed by Aradi et al 227 . DFTB parameter s for the Li - P - O - N chemical space were obtained using the recently develop ed TANGO (tight - binding approximation - enhanced global optimization) 228 method that allows fast and reliable parameterization by performing a small amount of DFT calculations. 6 . 2.1 Generation of c - Li 3 PO 4 and c - LiPON models - L i 3 PO 4 containing 128 atoms as initial c - Li 3 PO 4 structure. To build a c - LiPON cell, we replaced two O atoms by N atoms and removed two Li and two O atoms to balance the charge, so that the stoichiometry would be Li 2.875 PO 3.75 N 0.125 . Constant number of particles, pressure and temperature (NPT) ensemble was adopted to these two initi al cells at 300, 800, 1000, 1400 and 1500 K for 20 ps. Constant number, volume, and temperature (NVT) ensemble was implemented at 1400 K and 1500 K for 400 ps. 6 .2 .2 Generation of a - Li 3 PO 4 and a - LiPON models In order to obtain amorphous phases, we perform ed a heat - and - quench method as shown in Figure 6.1 . Initial structure guess was built by inserting 96 and 92 atoms in a cubic box by random packing, with stoichiome try Li 3 PO 4 and Li 2.83 PO 3.67 N 0.17 , respectively, and ~40% density 80 of the experimental referen ce for c - LiPON 79 , using the software Aten 229 . Then we annealed the cell at 1000 K for Li 3 PO 4 and 1400 K for LiPON , and quenched the st ructures at 300 K for equilibration, to obtain the amorphous phases. All t he heating and cooling steps were allow to equilibrate for 20 ps. NVT simulations were applied to a - Li 3 PO 4 and a - LiPON afterward for at least 200 ps , followed by equilibrated NPT at the certain temperature . Berendsen thermostat was used for NPT runs and N o se - Hoover thermostat was employed for NVT simulations. Time step was chosen to be 1 fs for both NPT and NVT ensembles. Figure 6 . 1 : Flowchart of generation of amorphous a - Li 3 PO 4 and a - LiPON . 81 6 . 3 Results and discussion 6 .3.1 Lattice parameters Table 6 . 1 : Calculated lattice parameters (a, b, c) of c - Li 3 PO 4 and c - LiPON at 300 K comparing to experiments. a(Å) b(Å) c(Å) c - Li 3 PO 4 (exp) 79 10.4612 6.1113 4.9208 c - L i 3 PO 4 (cal) 10.3267 6.0247 4.8499 Error (%) - 1.29 - 1.14 - 1.44 c - LiPON (exp) 79 10.4690 6.1153 4.9195 c - LiPON (cal) 10.3487 6.0376 4.8602 Error (%) - 1.15 - 1.27 - 1.21 We examined the lattice parameters obtained from NPT simulations at 300 K for c - Li 3 PO 4 and c - LiPON to validate our computational results (shown in Table 6. 1). Overa ll speaking the lattice parameters are in good agreement to experiments with about 1.1 - 1.5% underestimation. Volumes are also calculated for all four structures at different temperatures. Si nce we had different cell sizes for crystalline and amorphous phas es, for the sake of easy comparison, total volumes were normalized by dividing by 16 for crystalline phases and by 12 for amorphous phases. Linear relationship can be observed for crystallin e phase, while the amorphous phase had general increasing trend up on heating but kinks can be identified, e.g., amorphous Li 3 PO 4 at 800 K. Given that there were less atoms in amorphous LiPON cell than amorphous Li 3 PO 4 , the former had smaller volume than th e other. 6 .3.2 Exploration of LiPON structure As mentioned above, there is a controversy occurred in the field that whether the N t exist in LiPON . Spectroscopic works showed the evident that bridging N were presented in LiPON but it was difficult to distin guish whether they are N d or N t , or both. Here this computational work 82 enabled us to get atomic insight of LiPON structure. Fig ure 6. 2 showed an example of a - LiPON structure after NVT run at 1200 K, with N a and N d bonding environment highlighted. At all ot her temperatures we simulated, N a and N d coexisted in a - LiPON stru cture, and no N t was observed. - LiPON , two N d existed at 1500 K and two N a at 1400 K, which means the generation of N d need relatively large external energy. Figure 6 . 2 : Schema tic of a - LiPON structure at 1200 K from NVT run. To understand the role of N in a - LiPON structure, we also examined the bond distance (BD) and bond order (BO) of P - N and P - O. Several representative snapshots of trajectory from 1400 K NVT simulation of a - L iPON had been investigated. The Density Derived Electrostatic and Chemical (DDEC) scheme 151 c alculation was applied here to calculate the BO of P - N and P - O. In Fig ure 6. 3 we plotted BO as a function of BD of these pairs. BO of P - N a were larger than P - N d and P - O when the BD were the same, meaning that P - N a bonds were stronger than other two pairs. The longest BD of P - N a was 1.91 Å with BO = 1.09, while the shortest BD of this pair was 1.43 83 Å with BO = 2.17, which indicated that P - N a pair is dynamically changing b etween single bond and double bond. Same scenario can be observed for P - N d pair and P - O pair, but it is notable that the BD of P - N d can vary from 2.15 Å (BO = 0.53) to 1.46 Å (BO = 1.76). When N d was close to one P with BO close to 2, it was far from the P on the other side with BO close 1, indicating that two PO 3 radicals were linked by N d i n the form of P=N d - P, which was consistent with the hypothesis previously proposed by Wang et al 79 , although they claimed that only a small conc entration of this kind of units would contain in LiPON . This argument could be settled by increasing the population of N in future works. Figure 6 . 3 : Bond order as a function of bond distance of P - N, P - O and Li - X pairs. Similar analysis had been done for Li - X (X=N a , N d and O) pairs and was shown in Fig ure 6. 3 . Although Li atoms were likely to form weak ionic bond to other elements, here BO of Li - N a we re generally stronger than Li - O and Li - N d , suggesting that Li were less mobile when bonded to N a . In addition, we investigated the local environment of Li atoms in a - LiPON with the help of 84 electron localization function 230 (ELF) calculation. ELF maps of Li - X units with same BD of 2.12 Å were displayed in Fig ure 6.4 for comparison. Mushroom shapes electron sharing between Li - X were clear indicator of the existence of lone pair , while the electron shape was smaller and more l ocalized in Li - O than Li - N a and Li - N d . Less electron sharing of P - O than P - N a and P - N d can also be observed here. To further investigate the Li - X interactions, we selected 2000 snapshots from 200 ps trajecto r y at 1400 K and calculated the Li coordination e nvironment around these X atoms within cutoff distance of 2.3 Å, which was plotted in Fig ure 6. 5 . In the case of N a , more than half of the data located with 3 Li coordination, while 32.55% with 4 Li coordina t ion and 10.20% with two Li around. In contrast, the Li coordination around N d and O were more distributed, indicating that Li atoms were more mobile when bonded to N d and O. Figure 6 . 4 : ELF maps of (a) Li - N a , (b) Li - N d and (c) Li - O units (isosurface level of 0.82). 85 Figure 6 . 5 : Statistical results of Li coordination number around N a , N d and O. 6 .3.3 Self - diffusivity of different ato m ic groups Figure 6 . 6 : An example of as a function of Q 2 in Li diffusing at 1200 K for a - LiPON . Determination of Li self - diffusivity is critical for understanding mechanism of th e ionic 86 transport in Li - ion batteries 231 - 233 . Here we investigated the sel f - diffusivity of different atomic group by the incoherent density correlation I(Q,t) , as d iscussed in Section 3.3.2. We implemented this same method to all the atomic groups in Li 3 PO 4 and LiPON to examine the self - diffusivity and the diffusion mechanism. (a) (b) (c) (d) Figure 6 . 7 : of P and O for c - Li 3 PO 4 and c - LiPON at 1400 K and 1500 K. Figure 6. 6 is an example of as a function of Q 2 . Both Fickian model and SS model will be applied to extract the diffusivity D , residence time and mean square jump distance . The dynamics of atomic species P, O and N will be investigated first. A rapid decrease following by a 87 stable non - zero plateau indicated t he vibration motion in stead of diffusing. The of P and O with smallest Q value in c - Li 3 PO 4 and c - LiPON showed typical examples of vibrating and diffusing types, where we noticed that P and O were diffusing in the c - LiPON structure but remained vibrating in c - Li 3 PO 4 , a s shown in Fig ure 6. 7 . This implied that the addition of N atoms in c - LiPON helped the diffusion of P and O in the crystalline phase. The of N in c - LiPON was not included here but it showed typical di ffusion mode as well. Self - diffusivity of P, O and N were shown in Fig ure 6. 8 . Activation energy (E a ) was obtained by fitting the data to the Arrhenius law. Although the self - diffusivity of P, O and N can be calculated, the actual motion of these atoms could be correlated movement of a whole molecul e (e.g. PO 4 ), instead of individual jump. In c - LiPON , although P and N were diffusing according to their , their self - diffusivities were too low to get trustable estimation, and thus not included in the figure. It was notable that P and O were less mobile with higher activation energy in a - LiPON than in a - Li 3 PO 4 , especially P, which can be ascribed t o the existence of N d that bridged two P atoms together. The strong connection between N d and P due to their covalent bond, and thus a possible increase of network - strain energy 234 , remarkab ly reduced the mobility of P, as well as the O, which was bonded to P. In addition, in a - LiPON , P stopped diffusing below 1200 K, O and N stopped diffusing below 1000 K, suggesting the advantage of a - LiPON being a single ion (Li - only) conductor at low temp eratures. 88 Figure 6 . 8 : Self - diffusivity of P/O/N atomic groups from Li 3 PO 4 and L i PON as a function of temperature. Li self - diffusivity from different structures had been summarized in Fig ure 6. 9 , together with computational results from other work and experimental values. Our calculated Li self - diffusivity of a - Li 3 PO 4 was in good agre e experimental data measured by secondary ion mass spectroscopy 235 , with similar E a of 0.52 eV vs. 0.57 eV. Li self - diffusivity of a - Li 3 PO 4 obtained from Li et al. using neutral network (NN) potentials 236 was roughly one order of magnitude lower while the E a of 0.55 eV was c los e to ours. We also performed a DFT calculation using same a - L o PON structure at 1400 K employing Vienna Ab initio Simulation Package ( VASP ) 193 - 196 code for v alid ation. DFT result showed slightly higher Li self - diffusivity (1.24 x 10 - 4 cm 2 /s) than DFTB result (5.59 x 10 - 5 cm 2 /s) at 1400 K. In addition, a - LiPON was observed to have slightly lower Li self - diffusivity than a - Li 3 PO 4 . This is consistent with what we sta ted above, that the existence of N a in a - LiPON suppressed the Li diffusion due to stronger Li - N a bonds. c - LiPON had almost one order of magnitude higher Li self - diffusivity than c - Li 3 PO 4 (4.05 x 10 - 6 cm 2 /s vs. 6.24 x 10 - 7 cm 2 /s at 1400 K), while both h ad l arge E a (~2.70 eV). We also tried to extract Li diffusivity of crystalline phases below 1400 K and it turned out that Li barely moved. 89 The residence time and jump length of Li diffusion obtained from SS model fits were shown in Fig ure 6. 10 . Rapid incre ase of residence time can be observed for all structures and activation energy can be obtained by an Arrhenius fit. Activation energies extracted from residence time were comparable to the ones obtained from Li self - diffusivity. The jump lengths generally decr eased with increasing temperature, except the outlier for c - LiPON at 1500 K. Figure 6 . 9 : Li self - diffusivity of various structures as a function of temperature comparing to expe rime ntal results. (a) (b) Figure 6 . 10 : Residence time (a) and jump length (b) of Li diffusion. 90 6 .3.4 Ionic conductivity Ionic conductivity is known to be one of the most important properties for solid - state electrolyte materials. Here we extracted the ionic conductivity of Li 3 PO 4 and LiPON by coherent charge current density correlation function 237 - 238 , which includes the longitudinal and transverse components . The longitudinal component can be expressed as: , ( 38 ) while the transverse component as: , ( 39 ) where is the electron charge, the volume of the cell, the B oltzm ann constant, the temperature. is the collective charge current, which can be written as: , ( 40 ) where and is the charge and velocity of n t h ato m, respectively. The Fourier Transform (FT) then applied to the longitudinal and transverse component as: , ( 41 ) where is the frequency. The ionic conductivity is extr acted from the FT of transverse component as: . ( 42 ) Figure 6.11 reveals an example of the transverse component of coherent charge current density correlation function (a) and the c orresponding FT form (b) with three smallest Q values for a - LiPON at 1400 K . Unlike , which shows as a peak with zero intercept at zero frequency, has a non - zero intercept, meaning the ionic conductor 91 behavi or of a - Li PON . In the real case, we cannot access the zero frequency due to the limited cell size by applying Equation (42) . Thus, we average d the six outputs of smallest Q value from the smallest as the ionic conduct i vity. (a) (b) Figure 6 . 11 : (a) Time and (b) frequency domain for the transverse component of coherent charge current density function for three smallest Q at 1400 K for a - LiPON . Simulated resul ts were summarized in Fig ure 6. 1 2 , including several experimental data for comparison. The ionic conductivity of a - Li 3 PO 4 was comparable to the experimental results from Kuwata et al 235 and Yu et al 91 , except that the E a of 0.27 eV was slightly lower. Dyre et al. 239 energy, stating that the activat ion energy increased as temperature decreased, which might be accommodated in th i showed that a - LiPON had one order magnitude higher ionic conductivity than a - Li 3 PO 4 at lower temperature range, at higher temperatures (700 K to 1400 K), ionic conductivity of a - LiPON a nd a - Li 3 PO 4 were close according to our DFTB calculation results. This can be ascribed to the highly disorder of both structures, or the lack of charge carriers (Li + ), since a - LiPON had lower Li 92 concentration than a - Li 3 PO 4 . The excess Li concentration was proven to be effective by previous measurements and reached a maximum conductivity at 2.94 Li:P ratio 91 , 140 , 217 , 220 , 240 . Ionic conductivity of a - LiPON at 1400 K from DFT and DFTB showed good consistency, fur ther v alidating the reliability of our DFTB calculation. Experimental data from a similar composition Li 2.88 PO 3.73 N 0.14 conducted by Wang et al. 79 were plotted here for co mparison. c - LiPON revealed higher ionic conductivity and higher activation energy than c - Li 3 PO 4 , which was consistent with previous discussion in Fig ure 6 .7 , that only Li was diffusing in c - Li 3 PO 4 while all atomic groups were mobile in c - LiPON . Figure 6 . 12 : Calculated ionic conductivity of various structures as a function of temperature comparing to experimental data. 6 . 4 Summary DFTB is proven to be a powerful tool to deal w ith lar ge - scale MD simulation. In this work we firstly applied DFTB with new paramet ers for Li - P - O - N interactions to generate crystalline and amorphous LiPON structure and investigate the structural and dynamical properties of Li 3 PO 4 and LiPON . Both N a and N d can be observed in a - LiPON , but only in c - LiPON at higher 93 temperature, whereas n o evidence of N t existed. The bond order calculation indicated that N d incorporated with P in the form of P=N d - P, and Li interacted with N a more strongly than O or N d . In a ddition , self - diffusivity and ionic conductivity calculated by DFTB had been investi gated here. The addition of N atoms increased the mobility of P and O in crystalline phase. Li self - diffusivity of a - Li 3 PO 4 was in good agreement with experimental data, bu t less than one magnitude higher than results from other computational work. Jump length and residence time of Li diffusion from all four structures have been investigated. From our calculation, a - Li 3 PO 4 shared similar ionic conductivity with a - LiPON at th e high temperature range, which is likely due to the highly disordered structures for both. Further works need to be accomplished to verify some hypothesis raised by this work such as adjusting the N doping population and Li concentration in LiPON . The cap ab ility of DFTB demonstrated in this work also enable the feasibility of more complicated computational works, e.g., LiPON /Li metal interface study. 94 CHAPTER 7 : Preliminary study on LiPON /Li interphase 7 .1 Introduction The nobility of LiPON electrolyte comes from its ele ctrochemical stability contacting with Li metal, however, the origin remains i n conclusive . Sch wo bel et al. 2 41 pro posed the reaction at LiPON /Li interface with decomposition of LiPON and formation of Li 2 O, Li 3 N and Li 3 P, by X - ray photoemission spectroscopy (XPS). Another research by Zhu et al. 242 indicate d that the stability between electrolytes ( LiPON and Li 3 PS 4 ) and Li metal is due to the above side reaction products. DFTB is proven to be powerful dealing with larger unit cell, as well as modeling and predicting structural and dynamical properties of LiP ON in the previous chapter. Here in this chapter we will present some prel iminary study on LiPON /L i interphase using DFTB method , in order to monitor the reaction between LiPON /Li. 7 .2 Computational details Am orphous LiPON with composition Li 2.83 PO 3.67 N 0.17 obtained from Chapter 6 is applied as the bulk LiPON phase. Sur face ener gies of a, b and c direction of the a - LiPON with vacuum are calculated using the following equation by DFT: , ( 43 ) where is the total energy of the simulation box w ith the surface, the number of atoms and the energy per atom of bulk structure. Surface energies of a, b and c direction are obtained as 39.41 eV, 49.49 eV and 22.52 eV, respectively, indicating that c direction has the low est surface energy and thus be chosen as the contact surface with meta llic Li. 3x3x8 metallic Li slab is selected contacting with the LiPON phase. In order to accommodate the lattice parameters of 10.47 Å (a and b direction) of metallic Li slab, a stretche d 1x1x4 LiPON cell ( 10.47 95 Å on a, b direction and keep c direction unstretched) with vacuum on c direction is equilibrated at 450 K for 20 ps using NVT ensemble . The initial simulation super cell is composed by the equilibrated L i PON cell and a 3x3x8 metall ic Li (100) slab with 511 atoms in total, as shown in Figure 7.1. Metallic Li atoms are manually colored as dark green in order to distinguish against the Li atoms from LiPON , while color coding for other atoms are the same as Figure 6.1 . 2 Å vacuum layer is added between LiPON and metalli c Li on both sides. It is notable that a Li vacancy is created on purpose to mimic the realistic case. DFTB parameters for the Li - P - O - N chemical space were obtained using TANGO method. Then NVT ensemble is implemented to L iPON /Li supercell for 80 ps at 450 K to model the reaction and dynamics using cp2k 243 . Nose - Hoover thermostat is employed . Figure 7 . 1 : Initial supercell of LiPON /Li interface. 7.3 Results and discussion 7.3. 1 Average structure Average structure of the LiPON /Li supercell over 10 to 80 ps during the NVT run is shown in Figure 7.2 in order to evaluate the metallic Li layers. Additional simulation on metallic Li suggests that Li is bcc structure at 450 K under ou r computational setup. From Figure 7.2, i t is 96 notable that 3 layers of disordered metallic Li exist due to t he reaction with L i PON on both sides , while the middle layers are maintained as bcc structure . Figure 7 . 2 : Average structure over 10 to 80 ps during the NVT run. 7.3.2 Reaction at the interface Figure 7.3 depicts the MD simulation snapshots at 10 ps, 20 ps and 80 ps. No enormous change occurs during these period s , suggesting the reaction is accomplishe d with 10 ps , as well as the stability of the interface . It can be noted that the metallic Li atoms, which are close to interface, are attracted by the surface O atoms . No decomposition of L i PON , or formation of Li 2 O, Li 3 N and Li 3 P can be observed. In addi t ion, no Li ions exchange between these two phases can be identified during our simulation time. Sicolo et al. 224 proposed that rupture of P - N and P - O bonds occurred when using the model of amorphous Li 5 P 4 O 8 N 3 /Li metal, while in our case we cannot observe the rupture. 97 Figure 7 . 3 : Snapshots of LiPON /Li interphase at (a) 10 ps, (b) 20 ps and (c) 80 ps. 7.3.3 Li charge distribution In order to investigate the charge states of Li atoms in the model, we apply DDEC charge analys is 151 along the c axis, as shown in Figure 7.4. Since there is not much structural change during 10 to 80 ps from above, the snapshot structure at 20 ps is selected for the change analysi s a nd the following DOS calculation. The charge of Li atoms shows around 0.8 in the LiPON bulk phase and 0 in the met allic layer, which is consistent with our previous observation that the middle of metallic Li maintains its origin. While at the interface fro m 29 Å to 35 Å , the interface effect on the Li charge distribution can be identified , where Li charge differ from - 0.17 to 0.8, depending on the local environment and the distance from interface . At the interface region, the Li atoms 98 closer to interface sh ow higher charges. The charge distribution suggests that interface effect is confined into a 6 Å region between tw o phases. P (~1.3), O (~ - 0.9), N a (~ - 1.55) and N d (~ - 1.1) atoms do not show charge difference between the bulk and interface region. F igu re 7 . 4 : Variation of Li charge from DDEC charge calculation along the c axis. 99 7.3.4 Projected density of states Figure 7 . 5 : (a) A snapshot from MD at 20 ps with 3 selected layers (thickness of 6 Å) and two chosen Li atoms from metallic Li phase. (b) Projected density of states at selected layers and atoms. Fermi energy level is located at zero energy. The density of states projected at different layers and s e lected Li atoms and atoms are shown in Figure 7.5(b) by DFT . Li #1 and Li #2 are two metallic Li atoms with different distances from interface. Interface layer is chosen based on the previous charge analysis with thickness of 6 Å, and layer 1 and layer 2 a re two electrolyte layers with the same thickness, as shown in Figure 100 7.5(a). The Fermi level is located at zero energy. The metallic origin of Li #1 and Li #2 can be clearly identified regarding the states at Fermi energy level . It is notable that the pro jected DOS of interface also has states at the Fermi level, suggesting its metallic character . The metallic behavior of interface is likely due to the states induced by metallic Li. The electronic states at Fermi energy level re duce significantly beyond th e interface in layer 1 and layer 2, however, small amounts of states can be recorded, which is likely due to the limitation of this method as we take the whole supercell as unit cell and thus a uniform Fermi energy is obtained f or the whole cell. 101 CHAPTER 8 : F utu re work and c onclusions 8.1 Future work 8.1.1 Further investigation on t hermoelectric material tetrahedrite In this dissertation, we investigate d thermoelectric material tetrahedrite s on various aspects including their structures, low - temperature thermo ele ctric properties, origin of low lattice thermal conductivity and Cu movement in Cu - rich tetrahedrite. However, there are still a lot that remain s unknown for the tetrahedrite family w hich require further investigation. For instance, there are various po ssi ble dopants for Cu and Sb sites that have not been explored for tetrahedrite. Moreover, the origin of metal to semiconductor (MST) , accompanied with thermoelectric properties change, remains inconclusive. In addition, our study reveal ed the Cu movement in Cu - rich tetrahedrite , suggest ing its potential as an ionic conductor, however, the optimized composition and temperature is not well - understood yet. These future investigations might improve the performance of tetrahedrite and provide insight of materia ls design. 8.1. 2 Understanding the origin of LiPON /Li interphase Although the previous chapter provides some insights on the reaction of LiPON /LI inter phase , the origin of good electrochemical stability between the electrolyte - electrode remains unknown. D e tai led analysis on the structural and dynamic properties on the interface, such as Li ion migration and charge transfer across the interface , can help us un derstand the interface better. In addition, optimization of N composition is critical to the perform anc e of LiPON electrolyte. Experimental techniques such as surface morphology characterization, impedance spectroscopy are also helpful to extend our unders tandings and provide validation of our computational methods. Further study is needed on the reactio n a t the interface to provide research direction of improving battery applications. 102 8.1. 3 Exploration in cathode materials In addition to electrolyte materi als we discussed above , the selection and regulation of cathode materials are critical to the perfor man ce of Li - ion battery , especially in the application of electrical vehicles (EV). The number of EV has continued to increase rapidly around the world, which is estimated to be more than 125 million by the year of 2030 244 . Figure 1.7 shows some of the promi sing candidates for cathode materials . Among them, LiFePO 4 (LFP), Li(Ni 1 - x - y Mn x Co y )O 2 (NMC), and Li(Ni 1 - x - y Co x Al y )O 2 (NCA) are demonstrated to be outstanding and currently in commercial use. LFP (~170 Ah kg - 1 , ~3.45 V vs. Li/Li + ) has good electrochemical performance with low resistance, high current rating and long cycle life, while phosphate structure provides s tabi lit y against overcharging and heat, leads to a wide operating temperature range between - 245 . NMC (150 - 180 Ah kg - 1 depending on composition ) is another successful option due to its good performance by combing Ni and Mn, where Ni contributes to high specif ic en ergy and Mn offers the low internal resistance structure. Optimization of the percentage s of Ni, Mn and Co are attracting a lot of interests. NCA (~200 Ah kg - 1 , ~3.8 V vs. Li/Li + ) shares some similarities as NMC, but offers even higher specific energy , pow er densities and long life span 245 , however, safety concern exists. Improvement and regulation of current cathode materials, and exploration of new materials with advanced methods, are critical to further development of Li - ion batteries and EV industry. 8.2 Co nclus ions I n this thesis , structure and properties of two energy materials: tetrahedrite Cu 12 S 13 Sb 4 and L i PON have been investigated by experimental and computational techniques , such as DFT and DFTB, together with MD simulations. These advanced modeling techni ques a nd analysis methods employed in this thesis, are proven to be capable to explore the mechanism and predict materials 103 properties, which c an be i ntroduced to other material systems and provides insight on design of future energy material s. Low - temperat ure pr operties of Cu 12 Sb 4 S 13 attract many research interests due to the abrupt change of structural, electrical and thermoelectric properties at a specific metal - semiconductor transition (MST) temperature at around 85 K . In order to investigate the structu ral pr operties of Cu 12 Sb 4 S 13 , we use DFT - based molecular dynamics simulation along with neutron diffraction measurement . Cu 12 Sb 4 S 13 maintains as cubic structure according to our computational and experimental results. H owever, negative thermal expansion be low 50 K can be identified due to the interaction between Cu 12e and coordinated Sb atom . Low energy phonon vibrational mode can be observed. Thermoelectric properties such as electrical resistivity, Seebeck coefficient and electronic thermal conductivity a re als o examined below and above MST temperature comparing to experimental data. To understand the atomic dynamics of tetrahedrite, DFT - based first - principle molecular dynamics simulation is applied to investigate the N i and Zn co - doped tetrahedrite Cu 10.5 NiZn 0. 5 Sb 4 S 13 . The agreement between computed vibrational density of states and experimental inelastic neutron diffraction data at 300 K validate our simulation work. Anomalous phonon softening upon cooling can be explained by anharmonic rattling of Cu bet ween t wo Sb a toms inside the Sb[CuS 3 ]Cu atomic cage. The dynamic structure factors in the longitudinal and transverse direction obtained from coherent dynamics analysis suggest that the acoustic modes are confined in a small region of the scattering space, leadi ng to t he small thermal conductivity. Cu - rich tetrahedrite with composition Cu 1 4 Sb 4 S 13 is investigated in this thesis as well, focusing on the Cu transport properties. Structural instability can be speculated from the large fluctuation of lattice par ameter s durin g NPT run. Incoherent density correlation at 700 104 K indicates that part of Cu atoms in the structure keep vibrating, while the others are diffusing between different Cu sites. Transport properties of mobile Cu, such as self - diffusivity, res idence time and mean square jump lengt h can be extracted by applying Fickian model and SS model fit. Additionally, to visualize the actual motions of Cu atoms, nuclear density map is employed, which allows us to identify the Cu diffusion pathways. Migratio n energy barriers E m of different path ways E m for 24g 12d and 24g 12e diffusion paths are 1.32 eV and 1.56 eV, respectively, suggesting that when Cu atom locates at 24g site, it is more energetic favorable to jump to a 12d site t han to a 12e site. E ffect of simulation cell size and basis sets on the DFT - based molecular dyna mics simulation is examined using tetrahedrite Cu 10 Zn 2 Sb 4 S 13 thermoelectric as a model material. Three simulations , 1x1x1 cell with plane wave (PW) basis, 1x1x 1 cell with atomic orbital (AO) basis, and 2x2x2 cell with AO basis , are performed at 300 K to investigate the effect . Various structural and dynamic properties are extracted from three simulations, such as lattice parameters, partial and total vibrational density of states, phonon dispersion, heat capacity, and EXAFS spect ra . Comparison of these properties shows good agreement among these three simulations. We also investigate the 2x2x2 cell, which allows us to access the smaller Q range , where longitudina l and transverse sound/acoustic modes with linear dispersion can be i dentified. Sound velocities, bulk modulus and shear modulus can be extracted, showing good agreement with experimental data. Two low - energy quasi - localized vibrational modes can be observ ed from the coherent velocity correlation, which could scatter heat - c arrying acoustic modes by reducing their lifetime instead of reducing their speeds, leading to a small thermal diffusivity. With respect to LiPON electrolyte, we implemented advance DFTB method to investigate the structural and ionic conduction of Li 3 PO 4 and LiPON . New parameters for Li - P - O - N 105 interactions are applied to generate crystalline (c - ) and amorphous (a - ) Li 3 PO 4 and LiPON . Both apical N ( N a ) and doub ly - bridge N (N d ) can be observe d in a - LiPON , but only in c - LiPON at higher temperature, whereas no evidence of triply - bridge N ( N t ) exist s . Bond order calculation suggests that N d incorporated with P in the form of P=N d - P , and Li atoms have larger bond order with N a than with O or N d . M oreover, we examined the self - diffusivity of different atomic groups. With the addition of N, P and O atoms become more mobile in c - LiPON . The Li self - diffusivity of a - Li 3 PO 4 reveals good agreement with experimental data. Jump length and residence time o f Li diffusion from all four structures have been investigated as well. From the ionic conductivity calculation, we show that c - LiPON has one order magnitude higher ionic conductivity than c - Li 3 PO 4 , while a - Li 3 PO 4 shared similar ionic conductivity with a - LiP ON , which we ascribe to the highly disorder of both amorphous phases. This work provides insights of the amorphous structure of a - LiPON and the effect of N doping . In order to understand the origin of electroche mical stability between LiPON electrolyte and Li metal, we implement DFTB method to generate a LiPON /Li interface supercell and perform MD simulation on the supercell. From our simulation, only 3 layers of metallic Li involve the reaction with LiPON . Beyon d the 3 surface layers, the middle of metalli c Li atoms maintain as bcc structure based on our average structure analysis. Li are involving the reaction with LiPON on both sides, while the middle layers are maintained as bcc structure. The reaction at the interface occur within 10 ps, and n o decomposition of LiPON , or formation of Li 2 O, Li 3 N and Li 3 P can be observed. Li charge distribution is investigated, showing the thickness of interface is about 6 Å. The charges of P, O and N do not reveal differen ce between bulk and interface . Projected density of states of different layers and Li atoms are also 106 included. Metalli c behavior of the interface can be identified by the numbers of electronic states at Fermi energy level. 107 BIBL IO GRAPHY 108 BIBLIOGRAPHY 1. Petroleum, B., BP Statistical Review of World Energy Report. BP: Lon don, UK 2019 . 2. Energy, G., Global Energy & CO2 Status Report. International Energy Agency: Paris, France 2019 . 3. Change, I. C., The Physical Science Basis Cambridge, United Kingdom and New York, NY. USA: Cambridge University Press 2007, 2007 , 748 - 845. 4. Jenkinson, D. S.; Adams, D.; Wild, A., Model estimates of CO2 emissions from soil in response to global warming. Nature 1991, 351 (6324), 304. 5. Cherubini, F.; Peters , G. P.; Bernts en, T.; Strømman, A. H.; Hertwich, E., CO2 emissions from biomass combustion for bioenergy: atmospheric decay and contribution to global warming. Gcb Bioenergy 2011, 3 (5), 413 - 426. 6. Rosa, L. P.; Ribeiro, S. K., The present, past, and futu re contribution s to global warming of CO2 emissions from fuels. Climatic Change 2001, 48 (2 - 3), 289 - 307. 7. Dessler, A. E.; Parson, E. A., The science and politics of global climate change: A guide to the debate . Cambridge University Press: 2019. 8. Forman , C.; Muritala, I. K.; Pardemann, R.; Meyer, B., Estimating the global waste heat potential. Renewable and Sustainable Energy Reviews 2016, 57 , 1568 - 1579. 9. Elsheikh, M. H.; Shnawah, D. A.; Sabri, M. F. M.; Said, S. B. M.; Hassan, M. H.; Bashir, M. B. A.; Mohamad, M., A review on thermoelectric renewable energy: Principle parameters that affect their performance. Renewable and sustainable energy reviews 2014, 30 , 337 - 355. 10. Thackeray, M. M.; Wolverton, C.; Isaacs, E. D., Electrical energy storage for tra nsportation app roaching the limits of, and going beyond, lithium - ion batteries. Energy & Environmental Science 2012, 5 (7), 7854 - 7863. 11. Snyder, G. J.; Toberer, E. S., Complex thermoelectric materials. In materials for sustainable energy: a collection of peer - reviewed research and review articles from Nature Publishing Group , World Scientific: 2011; pp 101 - 110. 12. Chen, G.; Dresselhaus, M.; Dresselhaus, G.; Fleurial, J. - P.; Caillat, T., Recent developments in thermoelectric materials. International Materials Reviews 20 03, 48 (1), 45 - 66. 13. Tritt, T. M.; Subramanian, M., Thermoelectric materials, phenomena, and applications: a bird's eye view. MRS bulletin 2006, 31 (3), 188 - 198. 14. Alam, H.; Ramakrishna, S., A review on the enhancement of figure of merit from bulk to n ano - thermoelectric materials. Nano energy 2013, 2 (2), 190 - 212. 109 15. Price, S., The Peltier Effect and Thermoelectric Cooling. WWW - dokumentti. Saatavilla: http://ffden - 2 . phys. uaf. edu/212_s pring2007. web. dir/sedona_price/phys_212_webproj_peltier. html [v iitattu 20.5. 2017] 2007 . 16. Seebeck, T., Magnetische polarisation der metalle und erze durck temperatur - differenz, Abh. K. Akad. Wiss. Berlin 1823, 265 , 1823. 17. Riffat, S. B.; Ma, X., Th ermoelectrics: a review of present and potential applications. App lied thermal engineering 2003, 23 (8), 913 - 935. 18. Nolas, G.; Morelli, D.; Tritt, T. M., Skutterudites: A phonon - glass - electron crystal approach to advanced thermoelectric energy conversion applications. Annual Review of Materials Science 1999, 29 (1), 89 - 116. 19. Zhao, L. - D.; Lo, S. - H.; Zhang, Y.; Sun, H.; Tan, G.; Uher, C.; Wolverton, C.; Dravid, V. P.; Kanatzidis, M. G., Ultralow thermal conductivity and high thermoelectric figure of meri t in SnSe crystals. Nature 2014, 508 (7496), 373. 20. Heremans, J. P.; Jovovic, V.; Toberer, E. S.; Saramat, A.; Kurosaki, K.; Charoenphakdee, A.; Yamanaka, S.; Snyder, G. J., Enhancement of thermoelectric efficiency in PbTe by distortion of the electronic density of states. Science 2008, 321 (5888), 554 - 557. 21. Slack, G. A.; Rowe, D., CRC handbook of thermoelectrics. CRC press Boca Raton, FL: 1995. 22. Yang, J.; Stabler, F. R., Automotive Applications of Thermoelectric Materials. Journal of electronic mat erials 2009, 38 (7). 23. Matsubara, K. In Development of a high ef ficient thermoelectric stack for a waste exhaust heat recovery of vehicles , Twenty - First International Conference on Thermoelectrics, 2002. Proceedings ICT'02., IEEE: 2002; pp 418 - 423. 24. F leurial, J. - P., Thermoelectric power generation materials: Technol ogy and application opportunities. Jom 2009, 61 (4), 79 - 85. 25. Hendricks, T. J., Integrated thermoelectric thermal system resistance optimization to maximize power output in thermoelectric energy recovery systems. MRS Online Proceedings Library Archive 20 14, 1642 . 26. Tritt, T. M.; Böttner, H.; Chen, L., Thermoelectrics: Direct solar thermal energy conversion. MRS bulletin 2008, 33 (4), 366 - 368. 27. Min, G.; Rowe, D., Experimental evaluation of prototype thermoelectric domestic - refrigerators. Applied Energ y 2006, 83 (2), 133 - 152. 28. Vián, J.; Astrain, D., Development of a thermoelectric refrigerator with two - phase thermosyphons and capillary lift. Applied Thermal Engineering 2009, 29 (10), 1 935 - 1940. 110 29. Dai, Y.; Wang, R.; Ni, L., Experimental investigatio n and analysis on a thermoelectric refrigerator driven by solar cells. Solar energy materials and solar cells 2003, 77 (4), 377 - 391. 30. Phelan, P. E.; Chiriac, V. A.; Lee, T. - Y., Current an d future miniature refrigeration cooling technologies for high power microelectronics. IEEE Transactions on Components and Packaging Technologies 2002, 25 (3), 356 - 365. 31. Putra, N.; Sukyono, W.; Jo hansen, D.; Iskandar, F. N., The characterization of a ca scade thermoelectric cooler in a cryosurgery device. Cryogenics 2010, 50 (11 - 12), 759 - 764. 32. Choi, H. - S.; Yun, S.; Whang, K. - i., Development of a temperature - controlled car - seat system utilizing th ermoelectric device. Applied Thermal Engineering 2007, 27 (17 - 18), 2841 - 2849. 33. Riffat, S.; Qiu, G., Comparative investigation of thermoelectric air - conditioners versus vapour compression and absorption air - conditioners. Applied Thermal Engineering 2004, 24 (14 - 15), 1979 - 1993. 34. Lu, X.; Morelli, D. T.; Xia, Y.; Ozolins, V., Increasing the Thermoelectric Figure of Merit of Tetrahedrites by Co - Doping with Nickel and Zinc. Chemistry of Materials 2015, 27 (2), 408 - 413. 35. Lu, X.; Morelli, D. T.; Xia, Y.; Z hou, F.; Ozolins, V.; Chi, H.; Zhou, X.; Uher, C., High p erformance thermoelectricity in earth oÀ abund ant compounds based on natural mineral tetrahedrites. Advanced Energy Materials 2013, 3 (3), 342 - 348. 36. Heo, J.; Laurita, G.; Muir, S.; Subramanian, M. A. ; Keszler, D. A., Enhanced Thermoelectric Performance of Synthetic Tetrahedrites. Chemistry of Materials 2014, 26 (6), 2047 - 2051. 37. Momma, K.; Izumi, F., VESTA 3 for three - dimensional visualization of crystal, volumetric and morphology data. Journal of A pplied Crystallography 2011, 44 (6), 1272 - 1276. 38. Pfitz ner, A.; Evain, M.; Petricek, V., Cu12Sb4S13: A temperature - dependent structure investigation. Acta Crystallographica Section B: Structural Science 1997, 53 (3), 337 - 345. 39. Lai, W.; Wang, Y.; Morel li, D. T.; Lu, X., From bonding asymmetry to anharmonic r attling in Cu12Sb4S13 tetrahedrites: Whe n lone oÀ pair electrons are not so lonely. Advanced Functional Materials 2015, 25 (24), 3648 - 3657. 40. Suekuni, K.; Tomizawa, Y.; Ozaki, T.; Koyano, M., Systemat ic study of electronic and magnetic properties for Cu12 x TM x Sb4S13 (TM= Mn, Fe, Co, Ni, and Zn) tetrahedrite. Journal of Applied Physics 2014, 115 (14), 143702. 41. Bouyrie, Y.; Candolfi, C.; Ohorodniichuk, V.; Malaman, B.; Dauscher, A.; Tobola, J.; Len oir, B., Crystal structure, electronic band structure and high - temperature thermoelectric properties of Te - x 2.0). Journal of Materials Chemistry C 2015, 3 (40), 10476 - 10487. 111 42. Tablero, C., Electro nic and Optical Property Analysis of the Cu - Sb - S Tetrahed rites for High - Efficiency Absorption Devices. Journal of Physical Chemistry C 2014, 118 (28), 15122 - 15127. 43. Tippireddy, S.; Chetty, R.; Naik, M. H.; Jain, M.; Chattopadhyay, K.; Mallik, R. C., Ele ctronic and thermoelectric properties of transition metal substituted tetrahedrites. The Journal of Physical Chemistry C 2018, 122 (16), 8735 - 8749. 44. Suekuni, K.; Kim, F.; Takabatake, T., Tunable electronic properties and low thermal conductivity in synt 4, M= Ge, Sn). Journal of Applied Physics 2014, 116 (6), 063706. 45. Suekuni, K.; Tsuruta, K.; Kunii, M.; Nishiate, H.; Nishibori, E.; Maki, S.; Ohta, M.; Yamamoto, A.; Koyano, M., High - performance thermoelectric mi tetrahedrite. Journal of Applied Physics 2013, 113 (4), 043712. 46. Suekuni, K.; Tsuruta, K.; Ariga, T.; Koyano, M., Thermoelectric properties of mineral tetrahedrites Cu10Tr2Sb4S13 with low thermal conductivity. Applied Physics E xpress 2012, 5 (5), 051201. 47. Lu, X.; Morelli, D., The effect of Te substitution for Sb on thermoelectric properties of tetrahedrite. Journal of Electronic Materials 2014, 43 (6), 1983 - 1987. 48. Bouyrie, Y.; Candolfi, C.; Dauscher, A.; Malaman, B.; Lenoi r, B., Exsolution Process as a Route toward Extremely Low Thermal Conductivity in Cu12Sb4 p$ x Te x S13 Tetrahedrites. Chemistry of Materials 2015, 27 (24), 8354 - 8361. 49. Kumar, D. P.; Chetty, R.; Femi, O.; Chattopadhyay, K.; Malar, P.; Mallik, R., Thermoele ctric properties of Bi doped tetrahedrite. Journal of Electronic Materials 2017, 46 (5), 2616 - 2622. 50. Chetty, R.; DS, P. K.; Rogl, G.; Rogl, P.; Bauer, E.; Michor, H.; Suwas, S.; Puchegger, S.; Giester, G.; Mallik, R. C., Thermoelectric properties of a M n substituted synthetic tetrahedrite. Physical Chemistry Chemical Physics 2015, 17 (3), 1716 - 1727. 51. Weller, D. P.; Morelli, D. T., Rapid synthesis of zinc and nickel co - doped tetrahedrite thermoelectrics by reactive spark plasma sintering and mechanical alloying. J. Alloys Compd. 2017, 710 , 794 - 799. 52. Lu, X.; Morelli, D. T., Rapid synthesis of high - performance thermoelectric materials directly from natural mineral tetrahedrite. Mrs Communications 2013, 3 (3), 129 - 133. 53. Barbier, T.; Rollin oÀ Martinet, S.; Lemoine, P.; Gascoin, F.; Kaltzoglou, A.; Vaqueiro, P.; Powell, A. V.; Guilmeau, E., Thermoelectric Materials: A New Rapid Synthesis Process for Nontoxic and High oÀ Performance Tetrahedrite Compounds. Journal of the American Ceramic Society 2016, 99 (1), 51 - 56. 112 54. Meinel, A. B.; Meinel, M. P., Applied solar energy: an introduction. NASA STI/Recon Technical Report A 1977, 77 . 55. Boyle, G., Renewable energy. Renewable Energy, by Edited by Godfrey Boyle, pp. 456. Oxford University Press, May 2004. ISBN - 10: 0199261784. ISBN - 13: 9780199261789 2004 , 456. 56. Kenisarin, M.; Mahkamov, K., Solar energy storage using phase change materials. Renewable and sustainable energy reviews 2007, 11 (9), 1913 - 1965. 57. Zhao, H.; Wu, Q.; Hu, S.; Xu, H.; Rasmussen, C. N., Rev iew of energy storage system for wind power integration support. Applied energy 2015, 137 , 545 - 553. 58. Korpaas, M.; Holen, A. T.; Hildrum, R., Operation and sizing of energy storage for wind power plants in a market system. International Journal of Electr ical Power & Energy Systems 2003, 25 (8), 599 - 606. 59. Association, E. W. E., The economics of wind energy . EWEA: 2009. 60. Chum, H.; Faaij, A.; Moreira, J.; Junginger, H., Bioenergy. 2011 . 61. Bioenergy, I., Bioenergy a sustainable and reliable energy sou rce. International Energy Agency Bioenergy, Paris, France 2009 . 62. Rittmann, B. E., Opportunities f or renewable bioenergy using microorganisms. Biotechnology and bioengineering 2008, 100 (2), 203 - 212. 63. Goodenough, J. B.; Kim, Y., Challenges for recharg eable Li batteries. Chemistry of materials 2009, 22 (3), 587 - 603. 64. Zaghib, K.; Guerfi, A.; Hoving ton, P.; Vijh, A.; Trudeau, M.; Mauger, A.; Goodenough, J. B.; Julien, C., Review and analysis of nanostructured olivine - based lithium recheargeable batteri es: Status and trends. Journal of Power Sources 2013, 232 , 357 - 369. 65. Lu, Z.; Beaulieu, L.; Donabe rger, R.; Thomas, C.; Dahn, J., Synthesis, Structure, and Journal of The Electrochemical So ciety 2002, 149 (6), A778 - A791. 66. Bruce, P. G.; Freunberger, S. A.; Hardwick, L. J.; Tarascon, J. - M., Li O 2 and Li S batteries with high energy storage. Nature materials 2012, 11 (1), 19. 67. Xu, J.; Dou, S.; Liu, H.; Dai, L., Cathode materials for next generation lithium ion batteries. Nano Energy 2013, 2 (4), 439 - 442. 68. Taracson, J.; Armand, M., I ssues and challenges facing lithium ion batteries. nature 2001, 414 , 359 - 367. 69. Sawai, K.; Iwakoshi, Y.; Ohzuku, T., Carbon materials for lithium - ion (shu ttlecock) cells. Solid State Ionics 1994, 69 (3 - 4), 273 - 283. 113 70. Sandhya, C.; John, B.; Gouri, C., L ithium titanate as anode material for lithium - ion cells: a review. Ionics 2014, 20 (5), 601 - 620. 71. Goward, G.; Taylor, N.; Souza, D.; Nazar, L., The true crystal structure of Li17M4 (M= Ge, Sn, Pb) revised from Li22M5. Journal of alloys and compounds 200 1, 329 (1 - 2), 82 - 91. 72. Kamaya, N.; Homma, K.; Yamakawa, Y.; Hirayama, M.; Kanno, R.; Yonemura, M.; Kamiyama, T.; Kato, Y.; Hama, S.; Kawamoto, K., A lithi um superionic conductor. Nature materials 2011, 10 (9), 682. 73. Thangadurai, V.; Narayanan, S.; Pin zaru, D., Garnet - type solid - state fast Li ion conductors for Li batteries: critical review. Chemical Society Reviews 2014, 43 (13), 4714 - 4727. 74. Gwon, H.; Hong, J.; Kim, H.; Seo, D. - H.; Jeon, S.; Kang, K., Recent progress on flexible lithium rechargeable batteries. Energy & Environmental Science 2014, 7 (2), 538 - 551. 75. Zhang, S.; Ervin, M.; Xu, K.; Jow, T., Microporous poly (acrylonitrile - methyl methacryl ate) membrane as a separator of rechargeable lithium battery. Electrochimica Acta 2004, 49 (20), 333 9 - 3345. 76. QIU, W. - l.; YANG, Q. - h.; MA, X. - h.; FU, Y. - b.; ZONG, X. - f., Research on PEO - based dry solid polymer electrolytes for rechargeable lithium batter ies [J]. Chinese Journal of Power Sources 2004, 7 . 77. Morata - Orrantia, A.; García - Martín, S.; Morán , E.; Alario - Franco, M. Á., A New La2/3Li x Ti1 - x Al x O3 Solid Solution: Structure, Microstructure, and Li+ Conductivity. Chemistry of materials 2002, 14 ( 7), 2871 - 2875. 78. Bates, J.; Dudney, N.; Gruzalski, G.; Zuhr, R.; Choudhury, A.; Luck, C.; Robertso n, J., Electrical properties of amorphous lithium electrolyte thin films. Solid state ionics 1992, 53 , 647 - 654. 79. Wang, B.; Chakoumakos, B.; Sales, B.; Kw ak, B.; Bates, J., Synthesis, crystal structure, and ionic conductivity of a polycrystalline lithium - Li3PO4 structure. Journal of Solid State Chemistry 1995, 115 (2), 313 - 323. 80. Schnick, W.; Luecke, J., Synthesis and crystal structure of lithium phosphorus nitride Li7PN4: the first compound containing isolated PN4 - tetra hedra. Journal of Solid State Chemistry 1990, 87 (1), 101 - 106. 81. Daidouh, A.; Veiga, M.; Pico, C.; Martinez - Ripoll, M., A new polymorph of Li4P2O7. Acta Crystallographica Section C: Crystal Structure Communications 1997, 53 (2), 167 - 169. 82. Hermansen, C .; Mauro, J. C.; Yue, Y., A model for phosphate glass topology considering the modifying ion sub - network. The Journal of Chemical Physics 2014, 140 (15), 154501. 83. Du, Y. A.; Holzwarth, N., First - principles study of LiPON and related solid electrolytes. Physical Review B 2010, 81 (18), 184106. 114 84. Bunker, B. C.; Tallant, D. R.; Balfe, C. A.; Kirkpatrick, R. J.; Turner, G. L.; Reidmeyer, M. R., Structure of phosphorus oxynitride glasses. Journal of the American Ceramic Society 1987, 70 (9), 675 - 681. 85. Ma rchand, R.; L'Haridon, P.; Laurent, Y., Etude cristallochimique de LiPN2: Une structure derivée de la cristobalite. Journal of Solid State Chemistry 1982, 43 (2), 126 - 130. 86. Horstmann, S.; Irran, E.; Schnick, W., Phos phor (V) oÀ nitrid qfoÀ P3N5: Synthese ausg ehend von Tetraaminophosphoniumiodid und Kristallstrukturaufklärung mittels Synchrotron oÀ Pulver oÀ Röntgenbeugung. Zeitschrift f p› r anorganische und allgemeine Chemie 1998, 624 (4), 620 - 628. 87. West, W.; Whitacre, J.; Lim, J., Chemical stability enhancement of lithium conducting solid electrolyte plates using sputtered LiPON thin films. Journal of power sources 2004, 126 (1 - 2), 134 - 138. 88. Sagane, F.; Ikeda, K. - i.; Okita, K.; S ano, H.; Sakaebe, H.; Iriyama, Y., Effects of c urrent densities on the lithium plating morphology at a lithium phosphorus oxynitride glass electrolyte/copper thin film interface. Journal of Power Sources 2013, 233 , 34 - 42. 89. Albertus, P.; Babinec, S.; Lit zelman, S.; Newman, A., Status and challenges i n enabling the lithium metal electrode for high - energy and low - cost rechargeable batteries. Nature Energy 2018, 3 (1), 16. 90. Bates, J.; Lubben, D.; Dudney, N.; Hart, F., 5 Volt Plateau in LiMn2 O 4 Thin Film s. Journal of the electrochemical society 1995, 142 (9), L149 - L151. 91. Yu, X.; Bates, J.; Jellison, G.; Hart, F., A stable thin oÀ film lithium electrolyte: lithium phosphorus oxynitride. Journal of the electrochemical society 1997, 144 (2), 524 - 532. 92. Dud ney, N. J., Addition of a thin - film inorganic solid electrolyte (Lip on) as a protective film in lithium batteries with a liquid electrolyte. Journal of Power Sources 2000, 89 (2), 176 - 179. 93. Hohenberg, P.; Kohn, W., Inhomogeneous electron gas. Physical r eview 1964, 136 (3B), B864. 94. Kohn, W.; Sham, L. J., Self - consiste nt equations including exchange and correlation effects. Physical review 1965, 140 (4A), A1133. 95. Klenk, M. J.; Lai, W., Finite - size effects on the molecular dynamics simulation of fast - ion conductors: A case study of lithium garnet oxide Li7La3Zr2O12. S olid State Ionics 2016, 289 , 143 - 149. 96. Sarmiento - Perez, R.; Botti, S.; Marques, M. A., Optimized exchange and correlation semilocal functional for the calculation of energies of formati on. Journal of chemical theory and computation 2015, 11 (8), 3844 - 38 50. 97. Christopher, J. C., Essentials of Computational Chemistry. Theories and Models 2002 . 115 98. Bylander, D.; Kleinman, L., Good semiconductor band gaps with a modified local - density appr oximation. Physical Review B 1990, 41 (11), 7868. 99. Painter, G., I mproved correlation corrections to the local - spin - density approximation. Physical Review B 1981, 24 (8), 4264. 100. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized gradient approximat ion made simple. Physical review letters 1996, 77 (18), 3865. 101. S tephens, P. J.; Devlin, F.; Chabalowski, C.; Frisch, M. J., Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. The Journa l of physical chemistry 1994, 98 (45), 11623 - 11627. 102. Pribram - Jon es, A.; Gross, D. A.; Burke, K., Dft: A theory full of holes? Annual review of physical chemistry 2015, 66 , 283 - 304. 103. Zeng, Q.; Yu, A.; Lu, G., Multiscale modeling and simulation of po lymer nanocomposites. Progress in polymer science 2008, 33 (2), 191 - 269. 104. Berlin: 1991. 105. Fulde, P., Electron correlations in molecules and solids . Springe r Science & Business Media: 2012; Vol. 100. 106. Dai, J. In Investig ation of Li Ion Dynamics in Garnet Oxide Li6. 5La3Zr1. 5Ta0. 5O12 Using Molecular Dynamics Simulation , Meeting Abstracts, The Electrochemical Society: 2019; pp 174 - 174. 107. 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 Li 5 La 3 Ta 2 O 12 : A combined quasi - elastic neutron scattering and molecular dynamics study. Solid State Ionics 2017, 312 , 1 - 7. 1 08. Besenbacher, F.; Chorkendorff, I.; Clausen, B.; Hammer, B.; Mole nbroek, A.; Nørskov, J. K.; Stensgaard, I., Design of a surface alloy catalyst for steam reforming. Science 1998, 279 (5358), 1913 - 1915. 109. Kolmogorov, A.; Shah, S.; Margine, E.; Bialon, A.; Hammerschmidt, T.; Drautz, R., New superconducting and semicond ucting Fe - B compounds predicted with an ab initio evolutionary search. Physical review letters 2010, 105 (21), 217003. 110. Blouin, N.; Michaud, A.; Gendron, D.; Wakim, S.; Blair, E.; Neag u - Plesu, R.; Belletete, M.; Durocher, G.; Tao, Y.; Leclerc, M., Towa rd a rational design of poly (2, 7 - carbazole) derivatives for solar cells. Journal of the American Chemical Society 2008, 130 (2), 732 - 742. 111. Jain, A.; Shin, Y.; Persson, K. A., Computa tional predictions of energy materials using density functional theo ry. Nature Reviews Materials 2016, 1 (1), 15004. 116 112. Cohen, A. J.; Mori - Sánchez, P.; Yang, W., Challenges for density functional theory. Chemical reviews 2011, 112 (1), 289 - 320. 113. Elst ner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenh eim, T.; Suhai, S.; Seifert, G., Self - consistent - charge density - functional tight - binding method for simulations of complex materials properties. Physical Review B 1998, 58 (11), 7260. 114. Frauenheim, T.; Seifert, G.; Elstner, M.; Niehaus, T.; Köhler, C.; Amkreutz, M.; Sternberg, M.; Hajnal, Z.; Di Carlo, A.; Suhai, S., Atomistic simulations of complex materials: ground - state and excited - state properties. Journal of Physics: Condensed M atter 2002, 14 (11), 3015. 115. Gaus, M.; Cui, Q.; Elstner, M., DFTB3: e xtension of the self - consistent - charge density - functional tight - binding method (SCC - DFTB). Journal of chemical theory and computation 2011, 7 (4), 931 - 948. 116. Eschrig, H.; Bergert, I ., An optimized LCAO version for band structure calculations application to copper. physica status solidi (b) 1978, 90 (2), 621 - 628. 117. Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R., Construction of tight - binding - like potentials on t he basis of density - functional theory: Application to carbon. Physical R eview B 1995, 51 (19), 12947. 118. Foulkes, W. M. C.; Haydock, R., Tight - binding models and density - functional theory. Physical review B 1989, 39 (17), 12520. 119. Frauenheim, T.; Seif ert, G.; Elsterner, M.; Haj nal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R., A self oÀ consistent charge density oÀ functional based tight oÀ binding method for predictive materials simulations in physics, chemistry and biology. physica status solidi (b) 2000, 217 (1), 41 - 62. 120 . Seifert, G.; Porezag, D.; Frauenheim, T., Ca lculations of molecules, clusters, and solids with a simplified LCAO oÀ DFT oÀ LDA scheme. International journal of quantum chemistry 1996, 58 (2), 185 - 192. 121. Allen, M. P.; Tildesley, D. J., Computer simulation of liquids . Oxford university press: 2017. 122. LeSar, R., Introduction to computational materials science: fundamentals to applications . Cambridge University Press: 2013. 123. Tee, L. S.; Gotoh, S.; Stewart, W. E., Molecular parame ters for normal fluids. Lennard - Jones 12 - 6 Potential. Industrial & Engin eering Chemistry Fundamentals 1966, 5 (3), 356 - 363. 124. Wood, W.; Parker, F., Monte Carlo Equation of State of Molecules Interacting with the Lennard oÀ Jones Potential. I. A Supercritic al Isotherm at about Tw ice the Critical Temperature. The Journal of Chemical Physics 1957, 27 (3), 720 - 733. 117 125. Haftel, M. I., Surface reconstruction of platinum and gold and the embedded - atom model. Physical Review B 1993, 48 (4), 2611. 126. Stillinger, F. H.; Weber, T. A., Computer simulation of local order in condensed phases of silicon. Physical review B 1985, 31 (8), 5262 . 127. Carmesin, I.; Kremer, K., The bond fluctuation method: a new effective algorithm for the dynamics of polymers in all spatial dimensions. Macromolecules 1988, 21 (9), 2819 - 2823. 128. Van Gunsteren, W.; Berendsen, H., Algorithms for Brownian dynamics. Molecular Physics 1982, 45 (3), 637 - 647. 129. Hoogerbrugge, P.; Koelman, J., Simulating microscopic hydrodynamic phenomena with dis sipative particle dynamics. EPL (Europhysics Letters) 1992, 19 (3), 155. 130. Chen, S.; Doolen, G. D., Lattice Boltzmann met hod for fluid flows. Annual review of fluid mechanics 1998, 30 (1), 329 - 364. 131. Lee, B. P.; Douglas, J. F.; Glotzer, S. C., Filler - induced composition waves in phase - separating polymer blends. Physical Review E 1999, 60 (5), 5812. 132. Altevogt, P.; Ever s, O. A.; Fraaije, J. G.; Maurits, N. M.; van Vlimmeren, B. A., The MesoDyn project: software for mesoscale chemical engineering. Jo urnal of Molecular Structure: THEOCHEM 1999, 463 (1 - 2), 139 - 143. 133. Tucker III, C. L.; Liang, E., Stiffness predictions fo r unidirectional short - fiber composites: review and evaluation. Composites science and technology 1999, 59 (5), 655 - 671. 134. Burnet t, D. S., Finite element analysis: from concepts to applications. 1987 . 135. Guo, G.; Long, B.; Cheng, B.; Zhou, S.; Xu, P.; Cao, B., Three - dimensional thermal finite element modeling of lithium - ion battery in thermal abuse application. Journal of Power So urces 2010, 195 (8), 2393 - 2398. 136. Nasonova, D. I.; Verchenko, V. Y.; Tsirlin, A. A.; Shevelkov, A. V., Low - Temperature st ructure and thermoelectric properties of pristine synthetic tetrahedrite Cu12Sb4S13. Chemistry of Materials 2016, 28 (18), 6621 - 6627 . 137. Tanaka, H. I.; Suekuni, K.; Umeo, K.; Nagasaki, T.; Sato, H.; Kutluk, G.; Nishibori, E.; Kasai, H.; Takabatake, T., M etal Semiconductor Transition Concomitant with a Structural Transformation in Tetrahedrite Cu12Sb4S13. Journal of the Physical Socie ty of Japan 2015, 85 (1), 014703. 138. May, A. F.; Delaire, O.; Niedziela, J. L.; Lara - Curzio, E.; Susner, M. A.; Abernathy, D. L.; Kirkham, M.; McGuire, M. A., Structural phase transition and phonon instability in Cu 12 Sb 4 S 13. Physical Review B 2016, 93 (6), 064104. 118 139. Vaqueiro, P.; Gu lou, G.; Kaltzoglou, A.; Smith, R. I.; Barbier, T.; Guilmeau, E.; Powell, A. V., The influence of mobile copper ions on the glass - like thermal conductivity of copper - rich tetrahedrites. Chemistry of Materials 2017, 29 (9), 4080 - 4090. 140. Bates, J.; Dudney, N.; Gruzalski, G.; Zuhr, R.; Choudhury, A.; Luck, C.; Robertson, J., Fabrication and characterization of amorphous lithium electrolyte thin films and rechargeable thin - film batteries. Journal of power sources 1993, 43 (1 - 3), 103 - 11 0. 141. Di Benedetto, F.; Bernardini, G.; Cipriani, C.; Emiliani, C.; Gatteschi, D.; Romanelli, M., The distribution of Cu (II) and the magnetic properties of the synthetic analogue of tetrahedrite: Cu 12 Sb 4 S 13. Physics and chemistry of minerals 2005, 32 (3), 155 - 164. 142. Kresse, G.; Hafner, J., Ab initio molecular dynamics for liquid metals. Physical Review B 1993, 47 (1), 558. 143. Kresse, G.; Furthmüller, J., Efficiency of ab - initio total energy calculations for metals and semiconductors using a pla ne - wave basis set . Computational materials science 1996, 6 (1), 15 - 50. 144. Kresse, G.; Furthmüller, J., Efficient iterative schemes for ab initio total - energy calculations using a plane - wave basis set. Physical review B 1996, 54 (16), 11169. 145. Kresse, G.; Hafner, J., A b initio molecular - dynamics simulation of the liquid - metal amorphous - semiconductor transition in germanium. Physical Review B 1994, 49 (20), 14251. 146. Pizzi, G.; Volja, D.; Kozinsky, B.; Fornari, M.; Marzari, N., BoltzWann: A code for th e evaluation of t hermoelectric and electronic transport properties with a maximally - localized Wannier functions basis. Computer Physics Communications 2014, 185 (1), 422 - 429. 147. Pizzi, G.; Volja, D.; Kozinsky, B.; Fornari, M.; Marzari, N., An updated ver sion of BoltzWann : A code for the evaluation of thermoelectric and electronic transport properties with a maximally - localized Wannier functions basis. Computer Physics Communications 2014, 185 (8), 2311 - 2312. 148. Mostofi, A. A.; Yates, J. R.; Lee, Y. - S.; Souza, I.; Vander bilt, D.; Marzari, N., wannier90: A tool for obtaining maximally - localised Wannier functions. Computer physics communications 2008, 178 (9), 685 - 699. 149. Weller, D. P.; Stevens, D. L.; Kunkel, G. E.; Ochs, A. M.; Holder, C. F.; Morelli, D . T.; Anderson, M . E., Thermoelectric performance of tetrahedrite synthesized by a modified polyol process. Chemistry of Materials 2017, 29 (4), 1656 - 1664. 150. general fea tures. Zeitschrif t für Kristallographie - Crystalline Materials 2014, 229 (5), 345 - 352. 151. Manz, T. A.; Sholl, D. S., Improved atoms - in - molecule charge partitioning functional for simultaneously reproducing the electrostatic potential and chemical states i n periodic and no nperiodic materials. Journal of chemical theory and computation 2012, 8 (8), 2844 - 2867. 119 152. Lai, W.; Wang, Y.; Morelli, D. T.; Lu, X., From Bonding Asymmetry to Anharmonic Rattling in Cu12Sb4S13Tetrahedrites: When Lone - Pair Electrons Are Not So Lonely. Ad vanced Functional Materials 2015, 25 (24), 3648 - 3657. 153. Bouyrie, Y.; Candolfi, C.; Pailhes, S.; Koza, M.; Malaman, B.; Dauscher, A.; Tobola, J.; Boisron, O.; Saviot, L.; Lenoir, B., From crystal to glass - like thermal conductivity in cry stalline minerals . Physical Chemistry Chemical Physics 2015, 17 (30), 19751 - 19758. 154. Li, J.; Zhu, M.; Abernathy, D. L.; Ke, X.; Morelli, D. T.; Lai, W., First - principles studies of atomic dynamics in tetrahedrite thermoelectrics. APL Materials 2016, 4 ( 10), 104811. 155. Bullett, D.; Dawson, W., Bonding relationships in some ternary and quarternary phosphide and tetrahedrite structures:(Ag6M4P12) M6', Cu12+ xSb4S13 and Cu14 - xSb4S13, Ln6Ni6P17. Journal of Physics C: Solid State Physics 1986, 19 (29), 5837. 156. Bullett, D. W., Applications of atomic - orbital methods to the structure and properties of complex transition - metal compounds. Physics and Chemistry of Minerals 1987, 14 (6), 485 - 491. 157. Madsen, G. K.; Singh, D. J., BoltzTraP. A code for calculating band - structure d ependent quantities. Computer Physics Communications 2006, 175 (1), 67 - 71. 158. Ziman, J. M., Principles of the Theory of Solids . Cambridge university press: 1979. 159. Kresse, G.; Hafner, J., Ab initiomolecular dynamics for liquid metals. Physical Review B 1993, 47 (1), 558 - 561. 160. Kresse, G.; Hafner, J., Ab initiomolecular - dynamics simulation of the liquid - metal amorphous - semiconductor transition in germanium. Physical Review B 1994, 49 (20), 14251 - 14269. 161. Kresse, G.; Joubert, D., F rom ultrasoft pse udopotentials to the projector augmented - wave method. Physical Review B 1999, 59 (3), 1758. 162. Blöchl, P. E., Projector augmented - wave method. Physical review B 1994, 50 (24), 17953. 163. Fernandez - Alonso, F.; Price, D. L., Neutron Scatt ering . Academic P ress: 2013; Vol. 44. 164. Lu, X.; Morelli, D. T.; Xia, Y.; Zhou, F.; Ozolins, V.; Chi, H.; Zhou, X.; Uher, C., High Performance Thermoelectricity in Earth - Abundant Compounds Based on Natural Mineral Tetrahedrites. Advanced Energy Materials 2013, 3 (3), 342 - 348. 165. Baumert, J.; Gutt, C.; Shpakov, V.; Tse, J.; Krisch, M.; Müller, M.; Requardt, H.; Klug, D.; Janssen, S.; Press, W., Lattice dynamics of methane and xenon hydrate: Observation of symmetry - avoided crossing by experiment and theor y. Physical Revie w B 2003, 68 (17), 174301. 166. Tse, J.; Li, Z.; Uehara, K., Phonon band structures and resonant scattering in Na8Si46 and Cs8Sn44 clathrates. EPL (Europhysics Letters) 2001, 56 (2), 261. 120 167. Lara - Curzio, E.; May, A. F.; Delaire, O.; McGu ire, M. A.; Lu, X .; Liu, C. - Y.; Case, E. D.; Morelli, D. T., Low - temperature heat capacity and localized vibrational modes in natural and synthetic tetrahedrites. Journal of Applied Physics 2014, 115 (19), 193515. 168. Fan, X.; Case, E. D.; Lu, X.; Morelli , D. T., Room tem perature mechanical properties of natural - mineral - based thermoelectrics. J. Mater. Sci. 2013, 48 (21), 7540 - 7550. 169. Johnson, M. L.; Jeanloz, R., A Brillouin - zone model for compositional variation in tetrahedrite. American Mineralogist 1 983, 68 (1 - 2), 22 0 - 226. 170. Vaqueiro, P.; Guélou, G.; Kaltzoglou, A.; Smith, R. I.; Barbier, T.; Guilmeau, E.; Powell, A. V., The Influence of Mobile Copper Ions on the Glass - Like Thermal Conductivity of Copper - Rich Tetrahedrites. Chemistry of Materials 2 017, 29 (9), 4080 - 4090. 171. Makovicky, E.; Skinner, B. J., Studies of the sulfosalts of copper; VII, Crystal structures of the exsolution products Cu (sub 12.3) Sb 4 S 13 and Cu (sub 13.8) Sb 4 S 13 of unsubstituted synthetic tetrahedrite. The Canadian Mi neralogist 1979, 17 (3), 619 - 634. 172. Grimme, S.; Ehrlich, S.; Goerigk, L., Effect of the damping function in dispersion corrected density functional theory. Journal of computational chemistry 2011, 32 (7), 1456 - 1465. 173. Tran, F.; Stelzl, J.; Blaha, P., Rungs 1 to 4 of DFT Jac lattice constant, bulk modulus, and cohesive energy of solids. The Journal of chemical physics 2016, 144 (20), 204120. 174. Williams, G.; Watts, D. C., Non - symmetrical dielectric relaxation behavio ur arising from a simple empirical decay function. Transactions of the Faraday society 1970, 66 , 80 - 85. 175. Singwi, K.; Sjölander, A., Diffusive motions in water and cold neutron scattering. Physical Review 1960, 119 (3), 863. 176. Jobic, H.; Theodorou, D . N., Quasi - elastic neut ron scattering and molecular dynamics simulation as complementary techniques for studying diffusion in zeolites. Microporous and mesoporous materials 2007, 102 (1 - 3), 21 - 50. 177. Wang, Y.; Klenk, M.; Page, K.; Lai, W., Local Structu re and Dynamics of Lithi um Garnet Ionic Conductors: A Model Material Li5La3Ta2O12. Chemistry of Materials 2014, 26 (19), 5613 - 5624. 178. Jónsson, H.; Mills, G.; Jacobsen, K. W., Nudged elastic band method for finding minimum energy paths of transitions. 19 98 . 179. Henkelman, G.; Uberuaga, B. P.; Jónsson, H., A climbing image nudged elastic band method for finding saddle points and minimum energy paths. The Journal of chemical physics 2000, 113 (22), 9901 - 9904. 121 180. Henkelman, G.; Jónsson, H., Improved tange nt estimate in the nudge d elastic band method for finding minimum energy paths and saddle points. The Journal of chemical physics 2000, 113 (22), 9978 - 9985. 181. Palmer, D.; Conley, M., CrystalMaker. CrystalMaker Software, Bicester, England 2007 . 182. Kohn , W.; Sham, L. J., SELF - CONSISTENT EQUATIONS INCLUDING EXCHANGE AND CORRELATION EFFECTS. Phys. Rev. 1965, 140 (4A), 1133 - &. 183. Tran, F.; Stelzl, J.; Blaha, P., Rungs 1 to 4 of DFT Jacob's ladder: Extensive test on the lattice constant, bulk modulus, and cohesive energy of solid s. J. Chem. Phys. 2016, 144 (20). 184. Spiekermann, G.; Wilke, M.; Jahn, S., Structural and dynamical properties of supercritical H2O - SiO2 fluids studied by ab initio molecular dynamics. Chem. Geol. 2016, 426 , 85 - 94. 185. Miceli, G. ; Hutter, J.; Pasquarell o, A., Liquid Water through Density - Functional Molecular Dynamics: Plane - Wave vs Atomic - Orbital Basis Sets. J. Chem. Theory Comput. 2016, 12 (8), 3456 - 3462. 186. Ulian, G.; Tosoni, S.; Valdre, G., Comparison between Gaussian - type or bitals and plane wave ab initio density functional theory modeling of layer silicates: Talc Mg3Si4O10(OH)(2) as model system. J. Chem. Phys. 2013, 139 (20). 187. Suekuni, K.; Tsuruta, K.; Kunii, M.; Nishiate, H.; Nishibori, E.; Maki, S.; Ohta, M.; Yamamoto , A.; Koyano, M., High - p erformance thermoelectric mineral Cu12 - xNixSb4S13 tetrahedrite. J. Appl. Phys. 2013, 113 (4). 188. Lai, W.; Wang, Y.; Morelli, D. T.; Lu, X., From bonding asymmetry to anharmonic rattling: when lone - pair electrons are not so lonely. Adv. Funct. Mater. 2015 , 25 (24), 3648 - 3657. 189. Li, J. C.; Zhu, M. Z.; Abernathy, D. L.; Ke, X. L.; Morelli, D. T.; Lai, W., First - Principles Studies of Atomic Dynamics in Tetrahedrite Thermoelectrics. APL Materials 2016, 4 , 104811. 190. Zhou, F.; Niels on, W.; Xia, Y.; Ozolins , V., Lattice Anharmonicity and Thermal Conductivity from Compressive Sensing of First - Principles Calculations. Phys. Rev. Lett. 2014, 113 (18). 191. Suekuni, K.; Tomizawa, Y.; Ozaki, T.; Koyano, M., Systematic study of electronic a nd magnetic properties f or Cu12 - xTMxSb4S13 (TM=Mn, Fe, Co, Ni, and Zn) tetrahedrite. J. Appl. Phys. 2014, 115 (14). 192. Lu, X.; Morelli, D. T.; Wang, Y.; Lai, W.; Xia, Y.; Ozolins, V., Phase Stability, Crystal Structure, and Thermoelectric Properties of C u12Sb4S13 - xSex Sex Solid Solutions. Chem. Mater. 2016, 28 (6), 1781 - 1786. 193. Kresse, G.; Hafner, J., Abinitio molecular - dynamics for liquid - metals. Phys. Rev. B 1993, 47 (1), 558 - 561. 122 194. Kresse, G.; Furthmuller, J., Efficiency of ab - initio total energy calculations for metals and semiconductors using a plane - wave basis set. Comput. Mater. Sci. 1996, 6 (1), 15 - 50. 195. Kresse, G.; Furthmuller, J., Efficient iterative schemes for ab initio total - energy calculations using a plane - wave basis set. Phys. Rev. B 1996, 54 (16), 11169 - 11186. 196. Kresse, G.; Hafner, J., Ab - initio molecular - dynamics simulation of the liquid - metal amorphous - semiconductor transition in Germanium. Phys. Rev. B 1994, 49 (20), 14251 - 14269. 197. Kresse, G.; Joubert, D., From ultrasoft p seudopotentials to the p rojector augmented - wave method. Phys. Rev. B 1999, 59 (3), 1758 - 1775. 198. Blochl, P. E., Projector augmented - wave method. Phys. Rev. B 1994, 50 (24), 17953 - 17979. 199. VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chass aing, T.; Hutter, J., QUICKSTEP: Fast and accurate density functional calculations using a mix ed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167 (2), 103 - 128. 200. Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J., CP2K: atomist ic simulations of condensed matter systems. Wiley Interdisciplinary Reviews - Computational Mole cular Science 2014, 4 (1), 15 - 25. 201. Goedecker, S.; Teter, M.; Hutter, J., Separable dual - space Gaussian pseudopotentials. Phys. Rev. B 1996, 54 (3), 1703 - 1710. 202. Krack, M., Pseudopotentials for H to Kr optimized for gradient - corrected exchange - correl ation functionals. Theor. Chem. Acc. 2005, 114 (1 - 3), 145 - 152. 203. VandeVondele, J.; Hutter, J., Gaussian basis sets for accurate calculations on molecular syste ms in gas and condensed phases. J. Chem. Phys. 2007, 127 (11). 204. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77 (18), 3865 - 3868. 205. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalize d gradient approximation made simple (vol 77, pg 3865, 1996). Phys. Rev. Lett. 1997, 78 (7), 1 396 - 1396. 206. Martyna, G. J.; Tobias, D. J.; Klein, M. L., CONSTANT - PRESSURE MOLECULAR - DYNAMICS ALGORITHMS. J. Chem. Phys. 1994, 101 (5), 4177 - 4189. 207. Webb, S . M., SIXpack: a graphical user interface for XAS analysis using IFEFFIT. Phys. Scr. 2005, T11 5 , 1011 - 1014. 208. Ravel, B.; Newville, M., ATHENA, ARTEMIS, HEPHAESTUS: data analysis for X - ray absorption spectroscopy using IFEFFIT. Journal of Synchrotron Rad iation 2005, 12 , 537 - 541. 209. Johnson, N. E., X - Ray Powder Diffraction Data for Synthetic Var ieties of Tetrahedrite. Powder Diffr. 1991, 6 (1), 43 - 47. 123 210. Lara - Curzio, E.; May, A. F.; Delaire, O.; McGuire, M. A.; Lu, X.; Liu, C. - Y.; Case, E. D.; Morelli, D. T., Low - temperature heat capacity and localized vibrational modes in natural and synthetic tetrahedrites. J. Appl. Phys. 2014, 115 (19). 211. Rehr, J. J.; Kas, J. J.; Vila, F. D.; Prange, M. P.; Jorissen, K., Parameter - free calculations of X - ray spectr a with FEFF9. Phys. Chem. Chem. Phys. 2010, 12 (21), 5503 - 5513. 212. Hansen, J. P.; McDonald, I. R., Theory of Simple Liquids . Elsevier Ltd: 2013. 213. Shintani, H.; Tanaka, H., Universal link between the boson peak and transverse phonons in glass. Nat. Ma ter. 2008, 7 (11), 870 - 877. 214. Marchand, R., Nitrogen - containing phosphate glasses. Journal of Non - Crystalline Solids 1983, 56 (1 - 3), 173 - 178. 215. Day, D. E., Structural role of nitrogen in phosphate glasses. Journal of Non - Crystalline Solids 1989, 112 (1 - 3), 7 - 14. 216. Jacke, S.; Song, J.; Dimesso, L.; Brötz, J.; Becker, D.; Jaegermann, W., Tem perature dependent phosphorous oxynitride growth for all - solid - state batteries. Journal of Power Sources 2011, 196 (16), 6911 - 6914. 217. Fleutot, B.; Pecquenard, B.; Martinez, H.; Letellier, M.; Levasseur, A., Investigation of the local structure of LiPON thin films to better understand the role of nitrogen on their performance. Solid State Ionics 2011, 186 (1), 29 - 36. 218. Wang, B.; Kwak, B.; Sales, B.; Bates, J., Ionic conductivities and structure of lithium phosphorus oxynitride glasses. Journal of non - c rystalline solids 1995, 183 (3), 297 - 306. 219. Lacivita, V.; Westover, A. S.; Kercher, A.; Phillip, N. D.; Yang, G.; Veith, G.; Ceder, G.; Dudney, N. J., Resolving the amorphous structure of lithium phosphorus oxynitride (Lipon). Journal of the American Ch emical Society 2018, 140 (35), 11029 - 11038. 220. Lacivita, V.; Artrith, N.; Ceder, G., Structural and Compositional Factors That Control the Li - Ion Conductivity in LiPON Electrolytes. Chemistry of Materials 2018, 30 (20), 7077 - 7090. 221. Parr, R. G., Densi ty functional theory of atoms and molecules. In Horizons of Quantum Chemistry , Springer: 1980; pp 5 - 15. 222. Sicolo, S.; Albe, K., First - principles calculations on structure and properties of amorphous Li5P4O8N3 (LiPON). Journal of Power Sources 2016, 331 , 382 - 390. 223. Al - Qawasmeh, A.; Holzwarth, N., Li14P2O3N6 and Li7PN4: Computational study of two nitrogen rich cryst alline LiPON electrolyte materials. Journal of Power Sources 2017, 364 , 410 - 419. 124 224. Sicolo, S.; Fingerle, M.; Hausbrand, R.; Albe, K., Int erfacial instability of amorphous LiPON against lithium: A combined density functional theory and spectroscopic stud y. Journal of Power Sources 2017, 354 , 124 - 133. 225. Leung, K.; Pearse, A. J.; Talin, A. A.; Fuller, E. J.; Rubloff, G. W.; Modine, N. A., K inetics oÀ Controlled Degradation Reactions at Crystalline LiPON/LixCoO2 and Crystalline LiPON/Li oÀ M etal Interfaces. ChemSusChem 2018, 11 (12), 1956 - 1969. 226. Seifert, G., Tight - scheme. The Jour nal of Physical Chemistry A 2007, 111 (26), 5609 - 5613. 227. Aradi, B.; Hourahine, B.; Frauenheim, T., DFTB+, a sparse matrix - based implementation of th e DFTB method. The Journal of Physical Chemistry A 2007, 111 (26), 5678 - 5684. 228. Van den Bossche, M.; G r nbeck, H.; Hammer, B., Tight - binding approximation - enhanced global optimization. Journal of chemical theory and computation 2018, 14 (5), 2797 - 2807. 229. Cromie, S.; Elena, A.; Pashov, D.; Norman, S.; Forero - Martinez, N.; Smyth, M., Aten v1. 8 User Manu al Last Updated Tuesday, 03 December 2013. 230. Becke, A. D.; Edgecombe, K. E., A simple measure of electron localization in atomic and molecular syste ms. The Journal of chemical physics 1990, 92 (9), 5397 - 5403. 231. Coles, J.; Long, J., An ion - microprobe study of the self - diffusion of Li+ of lithium fluoride. Philosophical Magazine 1974, 29 (3), 457 - 471. 232. Takai, S.; Yoshioka, K.; Iikura, H.; Matsuba yashi, M.; Yao, T.; Esaka, T., Tracer diffusion coefficients of lithium ion in LiMn2O4 measured by neutro n radiography. Solid State Ionics 2014, 256 , 93 - 96. 233. Jambon, A.; Semet, M. P., Lithium diffusion in silicate glasses of albite, orthoclase, and obs idian composition: an ion - microprobe determination. Earth and Planetary Science Letters 1978, 37 (3), 445 - 450. 234. Anderson, O.; Stuart, D., Calculation of activation energy of ionic conductivity in silica glasses by classical methods. Journal of the Amer ican Ceramic Society 1954, 37 (12), 573 - 580. 235. Kuwata, N.; Lu, X.; Miyazaki, T.; Iwai, Y.; Tanabe, T.; Kawamura, J., Lithium diffusion coefficient in amorphous lithium phosphate thin films measured by secondary ion mass spectroscopy with isotope exchang e methods. Solid State Ionics 2016, 294 , 59 - 66. 236. Li, W.; Ando, Y.; Minamitani, E.; Watanabe, S., Stud y of Li atom diffusion in amorphous Li3PO4 with neural network potential. The Journal of chemical physics 2017, 147 (21), 214106. 237. Hansen, J. P.; M cDonald, I. R., Theory of simple liquids. Physics Today 1988, 41 , 89. 125 238. Dai, J.; Chen, Q.; Glossmann, T.; Lai, W., Comparison of interatomic potential models on the molecular dynamics simulation of fast - ion conductors: A case study of a Li garnet oxide Li7La3Zr2O12. Computational Materials Science 2019, 162 , 333 - 339. 239. Dyre, J. C.; Christensen, T.; Olse n, N. B., Elastic models for the non - Arrhenius viscosity of glass - forming liquids. Journal of non - crystalline solids 2006, 352 (42 - 49), 4635 - 4642. 240. Roh, N. - S.; Lee, S. - D.; Kwon, H. - S., Effects of deposition condition on the ionic conductivity and struc ture of amorphous lithium phosphorus oxynitrate thin film. Scripta materialia 1999, 42 (1). 241. Schwöbel, A.; Hausbrand, R.; Jaegermann, W., Interface reactions between LiPON and lithium studied by in - situ X - ray photoemission. Solid State Ionics 2015, 273 , 51 - 54. 242. Zhu, Y.; He, X.; Mo, Y., Origin of outstanding stability in the lithium solid electrolyte materials: insights from thermodynamic analyses based on first - principles calculations. ACS applied materials & interfaces 2015, 7 (42), 23685 - 23693. 243. Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J., cp2k: atomistic simulations of condensed matter systems. Wiley Interdisciplinary Reviews : Computational Molecular Science 2014, 4 (1), 15 - 25. 244. Miao, Y.; Hynan, P.; von Jouanne, A.; Yokoch i, A., Current Li - ion battery technologies in electric vehicles and opportunities for advancements. Energies 2019, 12 (6), 1074. 245. Hannan, M. A.; Hoqu e, M. M.; Hussain, A.; Yusof, Y.; Ker, P. J., State - of - the - art and energy management system of lithium - ion batteries in electric vehicle applications: Issues and recommendations. Ieee Access 2018, 6 , 19362 - 19378.