Thermal Conductivity of Solids from First-Principles Molecular

May 2, 2018 - Gotoh, Kunimitsu, Zhang, Lerner, Miyakubo, Ueda, Kim, Han, and Ishida. 0 (0),. Abstract: Diffusion of alkali metals in graphite layers i...
2 downloads 0 Views 1MB Size
Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

C: Energy Conversion and Storage; Energy and Charge Transport

Thermal Conductivity of Solids from FirstPrinciples Molecular-Dynamics Calculations John S Tse, Niall Joseph English, Ketao Yin, and Toshiaki Iitaka J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b00880 • Publication Date (Web): 02 May 2018 Downloaded from http://pubs.acs.org on May 5, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Thermal Conductivity of Solids from First-Principles Molecular-Dynamics Calculations

John S. Tse1*, Niall J. English2*, Ketao Yin1 and T. Iitaka3

1

Department of Physics and Engineering Physics University of Saskatchewan Saskatoon, Saskatchewan, Canada S7N 5E2

2

School of Chemical and Bioprocess Engineering University College Dublin, Belfield Dublin 4, Ireland 3

Computational Astrophysics Laboratory RIKEN 2-1 Hirosawa, Wako Saitama 351-0198, Japan

*

To whom correspondence should be addressed,

[email protected] and [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We describe the implementation of a computational scheme for estimating lattice thermal conductivities of ordered and disordered solids using the Einstein diffusion relationship and the energy moment, sampled from first-principles Born-Oppenheimer molecular dynamics employing Density Functional Theory.

The efficiency and accuracy of this procedure is

validated through calculations on of several compositionally diverse systems under different thermodynamic conditions, including MgO at elevated temperatures and pressures, as well as the potential efficient thermoelectric materials of doped Si46 and CoSb3 semiconductors. The results are found to be in good agreement with available experimental data covering a broad temperature and pressure range. Unlike some of current methods, no explicit calculation of high-order interplanar force constants or energy density is required. The method is most suitable for lowsymmetry crystalline, positionally-disordered and amorphous solids, particularly at high temperature.

ACS Paragon Plus Environment

Page 2 of 30

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction The thermal conductivity is an important thermophysical property of a material. Indeed, it is the fundamental quantity governing the performance of thermoelectric devices.1 The maximum efficiency for power generation and cooling is determined by the dimensionless figure of merit (ZT), which is proportional to the Seebeck coefficient and the electrical conductivity but inversely proportional to the thermal conductivity. Therefore, thermal conductivity is a critical parameter to be optimised. This subject is of topical interest, as thermoelectric generators have been reliable sources of electrical energy for remote terrestrial and extraterrestrial applications since the 1960s. More recently, there is increasing interest to harness waste heat, e.g. from industrial or geothermal sources or even vehicle exhausts, to provide clean electricity. Equally importantly, knowledge of relevant minerals’ thermal conductivity elevated pressure and temperature is sine qua non for understanding the rate of core heat loss and the Earth’s evolution.2

However, direct measurement of thermal conductivity under these extremes

conditions still poses formidable challenges. In recent years, there has been much efforts to develop computational methods for high-fidelity estimation of thermal conductivity from firstprinciples electronic-structure calculations. For systems that can be described by empirical pairwise interaction potentials, thermal conductivity is calculated most conveniently and accurately using equilibrium molecular dynamics (MD) employing the Green-Kubo relationship.3 This method requires the calculations of heat-flux correlation functions computed from the potential and kinetic energy of each particles and therefore, cannot be extended easily to first-principles electronic-structure-based MD (FPMD).3 A remedy is to employ interatomic potentials derived from fitting of the energies and forces from First-Principles electronic structure calculations. This approach is not entirely desirable as it is subjected to errors on the choice of the function representations and in the fitting procedure. In recent years, several methods have been developed to compute the thermal conductivity from First-Principles calculations.4 The most common computational method is to employ the phenomenological semi-classical Boltzmann transport equation (BTE).5 In this case, thermal conductivity is calculated from the kinetic theory mainly based on the phononrelaxation-time approximation,5–7 in which phonon scattering due to the anharmonic term of the potential expansion is calculated explicitly from the third- or higher-order interplanar force constants with the mode heat capacities, phonon frequencies, and group velocities determined

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 30

from harmonic lattice-dynamics calculations.8,9 The BTE method has been extended beyond the relaxation time approximation from solving the linearized BTE equation accurately.

10

Alternatively, finite-temperature anharmonic force constants can be extracted from analysis of intermediate scattering functions evaluated from an MD trajectory.11 Such a procedure is laborious and, indeed, prone to technical pitfalls. Moreover, since the temperature effect is only dependent on the phonon occupation of the Bose-Einstein distribution,8 it is not necessarily expected to provide a reliable result at temperatures close to melting or for intrinsically highly anharmonic systems, such as the recently discovered family of group-IV selenide thermoelectric materials.12 Moreover, explicit calculations of high-order force constants are computationally demanding13,14 and the current lattice-dynamics-based methods are only feasible for solids with high-symmetry space groups; applications to doped or disordered solids are very difficult, if not impossible. A non-equilibrium MD method in which the thermal conductivity is computed from the ratio of a heat flux to a temperature gradient has also been used in conjunction with FirstPrinciple MD.4,15 Recently, a formal expression for the adiabatic energy flux from DensityFunctional Theory (DFT) has been derived, which allows heat transport to be simulated using FPMD16 within the Green-Kubo formalism,17,18 and tested on liquid argon.16 A related method that avoids the “position-wrapping” problem has been also been proposed.19 In a previous study,20 employing a formulation based on the calculations of the energy momenta introduced by Kinaci et. al.,21 derived from pairwise additive many body potentials, we demonstrated that the calculation of the heat current in a solid can be circumvented and, instead, the thermal conductivity can be obtained directly from the time-dependent diffusion of energy momenta (i.e., thermal diffusion) of each of the system’s atoms, computed from first-principles MD via the Einstein-relationship as originally proposed by Helfand.

22,23

We found that such a

calculated thermal conductivity converges relatively rapidly with time and the size of the model system (albeit depending, in part, on the underlying phonons’ mean free path for sizeconvergence); in any event, this approach provides good agreement with results from the latticedynamics method. Motivated by this Einstein technique’s early encouraging results, in the present study, we apply this approach to study the thermal conductivity of several compositionally diverse systems under different thermodynamic conditions including crystalline and doped systems relevant to thermoelectric applications and a mineral under pressure. In more detail, to validate this method,

ACS Paragon Plus Environment

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

we investigate both temperature and pressure effects on the thermal conductivity of periclase (MgO) – considered very much a prototypical “test-bed” for thermal-conductivity calculations. Moreover, we gauge thermal conductivity of potential high-performance thermoelectric compounds, consisting of pure and Te-doped CoSb3 and Ba-doped silicon clathrate (Ba8Si46). As we shall present below, the results are in good accord with available experimental data. Theory The thermal conductivity of a material is defined via the Fourier law of heat conduction, in which the transported heat flux q is proportional to the thermal gradient (∂T/∂x),24 i.e.,

[1] where k is the thermal conductivity. The integrated heat current J is the rate of the energy transfer of each particle (εi) through a surface of cross-sectional area A (and via volume V) per unit time (vector normal to the surface) :3

[2] Defining the energy momentum R for the system as,21       

[3] we see that the heat current J is the just time derivative of R, i.e., J = dR/dt. Neglecting convection, the energy momentum for a particle k is21

    ∙  



[4] Eqn. (4) has been shown to be valid for N-body potentials.21 Although not derived from First Principles, this elegant construction for Rk in terms of the forces acting directly on the particle itself thus circumvents explicit calculation of the energy density.21 In addition, instead of using

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 30

the Green-Kubo formula to compute the thermal conductivity k from integrating the autocorrelation of the heat current J, this quantity can be computed using an equivalent Einstein relationship from the diffusion of the energy momentum, i.e., the mean squared displacement of the energy momentum, viz.20,21

κ

1 1 lim 〈 − 0    → 2

〉

[5] Since the instantaneous heat current J depends on the temporal positions, velocities and net forces of the particles, which are readily available in a classically-propagated MD calculation, it can be computed from MD using both empirical potentials and Density-Functional Theory (DFT) or other electronic-structure methods (i.e., first-principles Born-Oppenheimer MD). It is noteworthy that since the heat current is simply the time derivative of the energy moment, the advantages of the well-developed analytical “machinery” from Green-Kubo analysis of the heat current correlation functions, and its underlying intellectual infrastructure (e.g., for estimation of relaxation times characteristic of exponential-decay modes in thermal dissipation), are preserved fully, and to hand for the experienced practitioner to perform with skill. In the following. we will show that eqn. (4) is also applicable for FPMD calculations through selective examples on crystalline solids, doped and disordered crystals spanning a large thermal-conductivity range and under extreme temperature and pressure conditions. Computational details For DFT calculations on MgO, pure and Te-doped CoSb3 and Ba8Si46, the core electrons of the atoms were replaced by projected-augmented wave (PAW) potentials25 and a plane-wave expansion for the valence orbitals26 within the framework of the generalised-gradient approximation, 27

functional.

employing

the

Perdew-Burke–Ernzerhof

(PBE)

exchange-correlation

The Vienna ab initio simulation package (VASP) was used for all electronic and

first-principles Born-Oppenheimer molecular-dynamics calculations in the canonical (NVT) ensemble.28,29 For all cases, the integration time step was 2 fs. Additional calculations on MgO were performed using Vanderbilt ultra-soft pseudopotentials,30 and the results were compared with those obtained by PAW potentials. The supercell models used for the MD calculations were generated by replicating the respective primitive cells. For MgO,31 most calculations were

ACS Paragon Plus Environment

Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

performed on a 4×4×4 supercell. To test the convergence on the system size, a smaller 3×3×3 and a larger 5×5×5 supercell of MgO were also studied. For pristine and Te-doped CoSb3,32 and Ba8Si46 2×2×233 supercells were used resulting in 256 atoms for the former and 432 atoms for the latter. Integration in the Brillouin zone was sampled by the Γ point. Depending on the system and convergence, a trajectory consisted of 30,000 – 60,000 MD time steps (i.e. 60 – 120 ps) after thermal equilibration was required for thermal-conductivity calculations. We have taken special care to correct for periodic-image ‘jumps’ across the lattice imposed by periodic boundary conditions. In the post-simulation trajectory analysis on ‘wrapped’ coordinates, as an atom moves slightly outside as it vibrates about its equilibrium position without convection, it is allowed to stay outside naturally for the duration of that outside-box stay. In such a way, we employ ‘non-wrapped’ coordinates for the positions in the Einstein-approach analysis Periclase MgO Periclase is one of the most abundant mineral in the Earth lower mantle.34 In view of its importance, the thermal conductivity has been carefully determined over a broad temperature range under ambient pressure. Moreover, experimental data at high pressure, although not at the same quality of the ambient-pressure measurements, is also available. Therefore, this material is often chosen as a benchmark to validate new computational methods for thermal-conductivity prediction.4,8,11,35 Using the Einstein-MD procedure described above, we have computed the thermal conductivity of MgO at 100, 300 600, 900 and 1200 K under ambient pressure and at 300 K. The numerical results obtained using the 4×4×4 supercell with the ultra-soft pseudopotential for Mg an O atoms are presented in Table 1. Fig. 1 compares the theoretical thermal conductivity with experimental data and the results from previous anharmonic lattice-dynamics calculations.8 The calculated thermal conductivity of MgO at 300 K is 70.3±8.9 W/m/K; this is comparable to other theoretical predictions ranging from 65 – 111 W/m/K but is slightly higher than the consensus of reported experimental data, scattering from 36 to 70 W/m/K.36–40 The discrepancy may be attributed to imperfections of the density-functional treatment and the natural distribution of Mg and O isotopes in the experimentally-measured sample. The calculated thermal conductivity decreases from 153 W/m/K at 100K to 13.7 W/m/K at 1200K. The statistical errors were estimated from the averaged thermal conductivity computed from segments of the MD trajectory

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 30

at different time origins, but each of the same time duration. The errors become smaller at higher temperature as the sampling of the configuration space is more efficient. The MD results reproduced the expected trend for a crystalline material, due to the increase of anharmonicity underlying large-amplitude atomic motions, thus enhancing the phonon-phonon Umklapp scattering processes at high temperature. The absolute agreement between the computed thermal conductivity and experimental observations are very good over the entire temperature range. Both present and earlier lattice-dynamics calculations seem to overestimate the experimental values. The higher values are probably due to the neglect of the effect of thermal expansion in the theoretical calculations. If the volume increase at finite temperature was considered, the calculated thermal conductivity is expected to be lower, thus improving further the already-good agreement with experiment. Detailed comparison of the experimental thermal conductivities and those computed from anharmonic lattice-dynamics calculations shows the latter drop faster at increasing temperature. This is in contrast to the present MD results and experiment, in which the thermal conductivity remains fairly constant at high temperature. In the lattice dynamics approach, the harmonic phonon frequencies and the cubic anharmonic constants were computed at 0 K.8 These quantities were fixed and used in the calculations of thermal conductivity at finite temperature in which the temperature effect is derived from the variation of the phonon population governed by the Einstein distribution function. At high temperature, particularly near the melting point, higher-order force constants can, and do, become important and the harmonic frequencies need to be renormalised due to this anharmonic effect if quantitative accuracy is sought. In comparison, the phonon frequencies and the high-order effects are included implicitly in a MD calculation and therefore the Einstein analysis is more consistent with the experiments. To test system-size effects on the calculated thermal conductivities, we have repeated a calculation at 300 K using the larger 5×5×5 supercell model.

The computed thermal

conductivity was 73.3±9.2 W/m/K, which is within the error of the result from the smaller 4×4×4 model of 70.3±8.9 W/m/K. Therefore, the size dependence on this particular system is not too severe. These size-dependent results for MgO are also in general accord with smaller-system calculations for 3×3×3 supercells (66.8±8.6 W/m/K at 300 K), in that they are within estimated error-bar overlap. We have also studied the effect of using the PAW potential. The thermal conductivity of MgO at 300 K calculated from the PAW potentials using the 4×4×4 supercell is

ACS Paragon Plus Environment

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

73.6±8.0 W/m/K, in good agreement with the result obtained with ultra-soft potentials (vide supra). The measurement of thermal conductivity of a material under elevated pressure and temperature conditions is technically challenging. Most data were obtained from experimental measurements of thermal diffusivities and are converted to thermal conductivity from the heat capacity at constant pressure estimated from the equation of states and other elastic properties. Absolute measurements of thermal conductivity under extreme conditions have been made, but the results were subjected to be a variety of corrections due to the sample environment. The thermal conductivity of a MgO single crystal at room temperature has been measured to 60 GPa in a diamond-anvil cell40 by the transient thermo-reflectance spectroscopy,41 in which two subpicosecond pulses were generated from a single laser and used to pump (heat) and probe (reflectivity) the thermal diffusivity of a material. The calculated thermal conductivity of MgO at selected pressures at 300 K using the 4×4×4 supercell model are reported in Table 1. The thermal conductivities are found to increase with increasing pressure. At 137 GPa, the thermal conductivity of 357 ± 41 W/m/K is almost five times larger than the ambient-pressure value of 70.3 ± 8.9 W/m/K. To check the system-size dependence, an additional calculation using the larger 5×5×5 supercell was performed at 65 GPa. The result of 218 ± 31W/m/K is comparable to 211 ± 39 W/m/K computed from the smaller model. The calculated and measured thermal conductivities of MgO at 300 K as a function of pressure are compared in Fig. 2. In contrast to the temperature dependence, the statistical errors are larger at higher pressure, as the atoms’ motions are more restricted and longer MD runs may be necessary to achieve a more comprehensive sampling of the configuration space. However, the trend of increasing thermal conductivity with pressure is reproduced, although the absolute agreements are somewhat disappointing. Apparently, the thermal conductivities are overestimated somewhat by the theory and the deviation with experiment becomes larger with increasing pressure. However, one might expect an overestimation vis-à-vis experimental results due to the use in the present study of perfect, defect-free crystals, which facilitates phonons’ transfer of thermal energy. On the other hand, the present theoretical results are very similar to earlier anharmonic lattice-dynamics calculations.8 This discrepancy is puzzling, as there is no obvious reason that the MD method should perform less well under increasing pressure. A probable reason is that measured quantity in the experiment was the thermal diffusivity and then converted to thermal conductivity by

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

multiplying the estimated heat capacity corresponding to the pressure;40,41,42 however, the heat capacity at constant pressure is not known and cannot be measured directly. It is often estimated from the knowledge of the equation of state and other elastic properties from literature data.43,44 To test this speculation, we compare in Fig 3 the heat capacity at constant pressure (Cp) of MgO computed from the equation of state and quasi-harmonic (QHA) lattice dynamics45 at several pressures. The results show the Cp at ambient pressure are well reproduced by both methods, but the Cp estimated from the equation of state (i.e., so-termed “empirical” values) are lowered with respect to the anharmonic values at high pressures, and the deviation becomes larger at higher pressure. The ratios of the Cp at several pressures are tabulated in the inset of Fig. 3. If the “experimental” thermal conductivities were corrected by scaling with the appropriate (QHA/empirical) Cp ratios, the absolute magnitudes increased bringing the experimental results in better agreement with the theoretical Einstein-MD values. Ba-doped Si46 clathrate Thermal conductivity plays an important rôle in determining the thermoelectric performance of a material. The prevailing paradigm of “Phonon Glass Electron Crystal”46 for the design of advanced thermoelectric materials requires a potential candidate to possess good electrical conductivity and low thermal conductivity. One class of materials satisfy these criteria are metaldoped semi-conductor clathrates.47 It was anticipated that metal doping helps to increase the electron mobility and the localised ratting motions of the encaged metal atoms will enhance phonon-phonon scatterings and reduce the thermal conductivity.48 A representative material is silicon clathrates. Previously, crystalline silicon and silicon clathrates of types I (Si456) and II (Si136) structures were used to validate the Einstein-MD method for thermal-conductivity calculation presented here.20 From system-size dependence studies, it was found that the calculated thermal conductivities converged when the number of Si in the models is greater than ~300 atoms. Here, we studied the effect of Ba doping modelled by a 2×2×2 replicated supercell of type I silicon clathrate,33 consisting of 368 Si with both the small and large cavities filled completely by 32 Ba atoms (Ba8Si46) at ambient pressure. In view of the relatively large size of the model system and the use here of conventional O (N 3) DFT, MD calculations were performed only at 80 and 300 K. The results for the Ba-filled and empty Si46 reported previously are depicted in Fig. 4. The calculated lattice components to the thermal conductivities of Ba8Si46 at

ACS Paragon Plus Environment

Page 10 of 30

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

80 and 300 K are 1.8 ± 0.3 W/m/K and 3.3 ±0.5 W/m/K, respectively. The calculations reproduced correctly the glass-like behaviour of the thermal conductivity, with conductivity of the Ba-filled clathrate is lower than the empty clathrate and pure silicon. According to the resonant-scattering theory,48 in this temperature range, the thermal conductivity should be increasing slowly - corresponding to the “plateau” region of a disorder glass. A thermal conductivity of 1.7 ± 0.2 W/m/K has been reported from the measurement of a Ba8Si46 single crystal at 300 K by scanning thermal microscopy.49,50 The theoretical thermal conductivity reported here compares favourably with the experimental values. The thermal conductivity of Ba8Si46 has also been calculated from molecular dynamics employing interatomic potentials constructed by fitting the parameters to the trajectory of a DFT-based MD calculation.48 Using the Green-Kubo heat-flux time-correlation formalism, a thermal conductivity of 6 W/m/K was reported.48 The latter value is slightly higher than the present calculation, but within the same order of magnitude. The mechanism for phonon-phonon scatterings in filled clathrates has been explored thoroughly and will be elaborated here. In essence, the low-frequency Ba vibrations cross the acoustic branches of the Si clathrate. Due to the crystal symmetry, this crossing is symmetry-forbidden and resulted in strong interactions between the localised guest and dispersive lattice modes leading to exchange of heat energy between Ba and the Si lattice, thus reducing thermal conduction.51,52 Pure and Te-doped CoSb3 Skutterudite is another class of promising high-performance thermoelectric materials.53 Skutterudite is a naturally-occurring cobalt-arsenic mineral. Other variants with the arsenic replaced by other group-15 pnictogen elements, such as Te and Sb, can be synthesised. Unlike the clathrates, skutterudite does not have large open cavities. Instead, the Co atoms may be viewed as being encaged in the octahedra formed by the pnictogen elements. Single-crystal thermal conductivities of pristine CoSb3 from 100 – 700 K have been measured in two independent studies.54,55 A remarkably low thermal conductivity of 1.33–1.46 W/m/K at room temperature was reported. Furthermore, CoSb3 is a low-bandgap semiconductor, and a variety of elements can be doped substitutionally or into the interstitial sites to tune the electron concentration and mobility. Owing to these properties, skutterudites have generated substantial interest among all thermoelectric materials with applications at intermediate temperature.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 30

First-principles molecular-dynamics calculations in the NVT ensemble were performed at 150, 300, 600 and 800 K on a model constructed from 2×2×2 replication of the experimental unit cell. The respective computed thermal conductivities are 29.4 ± 5.5, 9.6 ± 1.7, 5.3 ± 0.8 and 5.3 ± 0.8 W/m/K. The theoretical and experimental data are compared in Fig. 5. The absolute agreement over a broad temperature range is exceedingly good.

More significantly, the thermal

conductivity is predicted to “plateau” around 5.3 W/m/K above 600 K. This is indeed observed by experiment. The thermal conductivity of crystalline CoSb3 has been calculated using abinitio lattice dynamics within the framework of the Boltzmann transport-equation method.54 The results, also shown in Fig. 4, show good agreement with experiment at low temperatures, but the discrepancy becomes larger when the temperature is above 650 K; the theoretical thermal conductivity continues to drop with increasing temperature, contrary to experimental observations. As described above (vide infra) a similar pattern was also observed in MgO. In addition, it is interesting to note that the thermal conductivities calculated from the local-density approximation (LDA) functional are in better agreement with experiment than those computed using the generalised-gradient approximation (GGA).56 Since the same GGA functional and the electronic-structure code were used in the present MD calculations, the difference is puzzling and further investigation is needed to resolve this problem. To study the effect of interstitial defects on the thermal conductivity, a Te atom was inserted in the empty tetrahedral site in the supercell model. This system has a stoichiometry of Te0.16,CoSb3. Since the doped system is metallic, only the lattice contribution to the thermal conductivity is calculated. At 300 K, a lower thermal conductivity of 4.05 ± 0.05 W/m/K than pure CoSb3 of 9.6 ± 1.66 W/m/K was predicted. Therefore, similar to Ba-doped Si clathrate (vide supra), the added Te atom serves as a “rattler” in scattering the heat-carrying lattice phonons and reducing the thermal conduction. However, the calculated thermal conductivity computed at 80 K of 17.2 ± 3.8 W/m/K is higher than at 300 K showing the Te-doped CoSb3 still retains crystalline thermal behavior. To investigate the substitutional effect, randomly-chosen Sb atoms in the 2×2×2 supercell of CoSb3 were replaced by Te to generate a model with stoichiometry Co4 Sb12-xTex (x = 0.2). The calculated thermal conductivities at 300 and 600 K are 4.62 ± 0.56 W/m/K and 3.65 ± 0.41 W/m/K, respectively. Experimental thermal-conductivity data on well-characterised high-quality single-crystal Te-doped CoSb3 samples are not available.

ACS Paragon Plus Environment

Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Moreover, since the atomic radius of Te (2.50 Å) is similar to that of Sb (2.45 Å), and they also have a very similar electronegativity (2.1 and 2.05 for Te and Sb, respectively), it is often difficult to ascertain the positions of Te atoms whether they are located substitutionally in the Sb sites or in the interstitial sites in the doped samples. The theoretical results are compared to the measured values on pressed power samples with the stiochiometry x = 0.1 and 0.3 in Fig. 6.57 It is gratifying that the predicted “crystalline-like” temperature trend and magnitudes of the thermal conductivities are comparable to the measured values considering the potential uncertainties with the experimental data. Conclusions In this paper, a method on computing the thermal conductivity of a solid using first-principles MD from the energy momenta20,21 is validated from consistently good numerical agreements with experiments on a variety of crystalline systems. The method reproduced the thermal conductivity of periclase (MgO) over a large temperature and pressure range. At high pressure, the thermal conductivity of a sample cannot be determined directly but this information can be converted from the measurement of the thermal diffusivity provided the heat capacity at constant pressure (Cp) is known. It is shown that the estimated Cp based on the equation of state is not reliable at high pressure. Better agreement between theoretical and “experimental” thermal conductivity may be achieved if the theoretical heat capacity calculated from quasi-harmonic approximation was used in the conversion. We also obtained good absolute agreement with experiments for the thermal conductivities of crystalline Ba8Si46 and CoSb3, two promising highperformance thermoelectric materials. In particular, the glass-like behaviour observed in Ba8Si46 is reproduced correctly. For pure and Te-doped CoSb3, the thermal conductivities are predicted to be relatively low; the crystalline behaviour is maintained. The computational algorithm presented here is valid only for solids: the convection contribution to thermal conduction is neglected. The potential problem associated with the latter has been resolved recently using a Berry-phase approach.16 When a crystal is at thermodynamic equilibrium (in this case, canonical ensemble, NVT), only the potential term contributes to the energy momenta, as the convection term is negligible. In principle, one cannot partition the energy momenta into individual atomic component. However, taking the equilibrium potential energy of the system as the reference, the temporal changes in the total energy can be

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

approximated as a sum of the perturbation on the individual atoms. Following this assumption, the energy change on each particle is simply the displacement times the net force acting on that atom, which is also physically intuitive. The forces can be computed via the Hellman-Feynman theorem, as were done here. The rate of change of the energy momenta can then be approximate as a sum of the individual atoms as in eqn. (4) of the manuscript. Therefore, eqn. (4) is a good first approximation to the full quantum mechanical treatment. Nevertheless, the computational scheme described here is accurate and simple to interface with existing first-principles electronic-structure MD codes. Comparing to current and widely-used methods based on the Boltzmann transport equation and relaxation times computed from anharmonic force constants or derived from MD trajectories through calculation of the intermediate scattering functions, the present method is more suitable and efficient to study both low-symmetry crystalline, positionally-disordered and amorphous solids. In the conventional BTE method based on the computation of the calculations of the cubic (or higher force constants), numerical finite differences of the forces computed on atomic positions displaced from equilibrium are needed. The displacements were selected from the symmetry of the crystal and usually only encompassed the nearest neighbour atoms. Therefore, for high symmetry crystals, such as MgO, the number of displacements are small. For non-crystalline solids modeled by supercells, since there is no crystal symmetry to reduce the number of atomic displacements, the number of calculation needed are very large and becoming unmanageable. Moreover, the number of displacements increases exponentially with the cutoff distance. The MD method presented here eliminates ipso facto this force-constant complexity. Moreover, the MD method can be applied to high temperature where high-order anharmonic effects certainly cannot be neglected if accuracy is desired. Furthermore, heat fluxes can be computed from the numerical derivative of energy momenta, permitting analysis of heat-flow mechanisms via fitting to a Fourier series of relaxation times.52 Alternatively, the contribution from individual phonon modes can be analysed from the MD trajectory by projecting to the energy-density correlation or intermediate scattering functions. It should be noted that quantum corrections to the calculated classical thermal conductivity wasere not implemented here. This issue is well known but seldom discussed amongst the practitioners of thermal -conductivity calculations. The root of the problem is due to the phonon -occupation number following Bose–Einstein statistics at finite temperature: as below the Debye

ACS Paragon Plus Environment

Page 14 of 30

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

temperature, not all the phonon modes are as much populated as in Maxwell–Boltzmann statistics. In a study on the thermal conductivity of pure and doped diamond from MD using empirical force-field and the thermal conductivities evaluated from the Green-Kubo heat flux relationship,58 it was found that the thermal conduction calculated from the classical correlation functions are similar to that computed with quantum corrections. Apparently, the quantum effect is not so important in this highly harmonic system. A detailed analysis of the contribution of the quantum corrections was performed on silicon59 using the phonon -group velocities and relaxation times calculated from lattice dynamics and thermal conductivity computed from the Boltzmann transport equation under the relaxation-time approximation. In this work, it is shown that the classical and quantum corrected thermal conductivities are in good agreement above the Debye temperature (674 K), but the classical results are lower at lower temperatures. Nevertheless, acceptable agreements were found at temperature as low as 200 K, about one-third of the Debye temperature. Therefore, near and above the Debye temperature, the classical approximation is acceptable. In the examples chosen for this study, the Debye temperature for MgO, Ba8Si46 and CoSb3 are, respectively 743 K,60 370 K61 and 306 K62 Except for MgO, the Debye temperatures of Ba8Si46, CoSb3 and doped Co3Sb12 are close or lower than 300 K, the temperature chosen for the MD simulations. If we can infer from the study on silicon60, for MgO at 300K, which is about 0.4 of the Debye temperature, the classical results maybe still be acceptable. On the other hand, the calculated thermal conductivity of BaSi46 at 80 K, which is about 0.2 of the Debye temperature, may be slightly underestimated; however, but the qualitative trend of almost constant thermal conductivity from 80 -300 K is correctly predicted. With current institutional high-performance computing capability, it is feasible to perform firstprinciples MD on several hundred atoms for durations of up to 100-200 ps. The method introduced here will probably become routine in the near future due to anticipated advances in more efficient algorithms for (linear-scaling) electronic-structure calculations and increasing computing power. Indeed, in our initial study on Si clathrates, we have already performed FPMD on model systems with more than 1000 atoms employing the O (N) method.20,63

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Acknowledgments JST wished to thank NSERC Canada for a discovery grant. Calculations were performed at the WestGrid Computing Facilities supported by ComputeCanada and on the HOKUSAI system at RIKEN (Japan) partly supported by MEXT as “Exploratory Challenge on Post-K computer” (Frontiers of Basic Science: Challenging the Limits).

ACS Paragon Plus Environment

Page 16 of 30

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

References (1)

Snyder, G. J.; Toberer, E. S. Complex Thermoelectric Materials. Nat. Mater. 2008, 7, 105–114.

(2)

A.M. Hofmeister, J. M. B. and M. P. Treatise in Geophysics; Price, G. D., Ed.; 2007.

(3)

Ladd, A. J. C.; Moran, B.; Hoover, W. G. Lattice Thermal Conductivity: A Comparison of Molecular Dynamics and Anharmonic Lattice Dynamics. Phys. Rev. B 1986, 34, 5058.

(4)

Stackhouse, S.; Stixrude, L. Theoretical Methods for Calculating the Lattice Thermal Conductivity of Minerals. Rev. Mineral. Geochemistry 2010, 71, 253–269.

(5)

McGaughey, A. J. H.; Kaviany, M. Quantitative Validation of the Boltzmann Transport Equation Phonon Thermal Conductivity Model under the Single-Mode Relaxation Time Approximation. Phys. Rev. B 2004, 69, 94303.

(6)

Callaway, J. Model for Lattice Thermal Conductivity at Low Temperatures. Phys. Rev. 1959, 113, 1046.

(7)

Narayanamurti, V.; Pohl, R. O. Tunneling States of Defects in Solids. Rev. Mod. Phys. 1970, 42, 201.

(8)

Tang, X.; Dong, J. Lattice Thermal Conductivity of MgO at Conditions of Earth’s Interior. Proc. Natl. Acad. Sci. 2010, 107, 4539–4543.

(9)

Togo, A.; Chaput, L.; Tanaka, I. Distributions of Phonon Lifetimes in Brillouin Zones. Phys. Rev. B 2015, 91, 94306.

(10)

Li, W.; Mingo, N. Lattice Dynamics and Thermal Conductivity of Skutterudites CoSb3 and IrSb3 from First Principles: Why IrSb3 Is a Better Thermal Conductor than CoSb3. Phys. Rev. B 2014, 90, 94302.

(11)

De Koker, N. Thermal Conductivity of MgO Periclase from Equilibrium First Principles Molecular Dynamics. Phys. Rev. Lett. 2009, 103, 1–4.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(12)

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 Merit in SnSe Crystals. Nature 2014, 508 (7496), 373–377.

(13)

Tadano, T.; Gohda, Y.; Tsuneyuki, S. Anharmonic Force Constants Extracted from FirstPrinciples Molecular Dynamics: Applications to Heat Transfer Simulations. J. Phys. Condens. Matter 2014, 26, 225402.

(14)

McGaughey, A.; Larkin, J. M. Predicting Phonon Properties from Equilibrium Molecular Dynamics Simulations. Ann. Rev. Heat Transf. 2014, 17, 49–87.

(15)

Müller-Plathe, F. A Simple Nonequilibrium Molecular Dynamics Method for Calculating the Thermal Conductivity. J. Chem. Phys. 1997, 106, 6082–6085.

(16)

Marcolongo, A.; Umari, P.; Baroni, S. Microscopic Theory and Quantum Simulation of Atomic Heat Transport. Nat. Phys. 2016, 12, 80–84.

(17)

Green, M. S. Markoff Random Processes and the Statistical Mechanics of Time‐ dependent Phenomena. II. Irreversible Processes in Fluids. J. Chem. Phys. 1954, 22, 398– 413.

(18)

Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. J. Phys. Soc. Japan 1957, 12, 570–586.

(19)

Kang, J.; Wang, L.-W. First-Principles Green-Kubo Method for Thermal Conductivity Calculations. Phys. Rev. B 2017, 96, 20302.

(20)

English, N. J.; John, S. T. Equilibrium Born-Oppenheimer Molecular-Dynamics Exploration of the Lattice Thermal Conductivity of Silicon Clathrates. Comput. Mater. Sci. 2017, 126, 1–6.

(21)

Kinaci, A.; Haskins, J. B.; Çağın, T. On Calculation of Thermal Conductivity from Einstein Relation in Equilibrium Molecular Dynamics. J. Chem. Phys. 2012, 137, 14106.

(22)

Helfand, E. Transport Coefficients from Dissipation in a Canonical Ensemble. Phys. Rev. Lett. 1960, 119, 1–9.

ACS Paragon Plus Environment

Page 18 of 30

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(23)

Zwanzig, R. Time-Correlation Functions and Transport Coefficients in Statistical Mechanics. Annu. Rev. Phys. Chem. 1965, 16, 67–102.

(24)

McQuarrie, D. A. Statistical Mechanics, Harper &. Row Publishers 1976.

(25)

Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector AugmentedWave Method. Phys. Rev. B 1999, 59, 1758.

(26)

Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, T. A.; Joannopoulos, J. D. Iterative Minimization Techniques for Ab Initio Total-Energy Calculations: Molecular Dynamics and Conjugate Gradients. Rev. Mod. Phys. 1992, 64, 1045.

(27)

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865.

(28)

Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186.

(29)

Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Open-Shell Transition Metals. Phys. Rev. B 1993, 48, 13115.

(30)

Vanderbilt, D. Soft Self-Consistent Pseudopotentials in a Generalized Eigenvalue Formalism. Phys. Rev. B 1990, 41, 7892.

(31)

Hazen, R. M. Effects of Temperature and Pressure on the Cell Dimension and X-Ray Temperature Factors of Periclase. Am. Mineral. 1976, 61, 266–271.

(32)

Schmidt, T.; Kliche, G.; Lutz, H. D. Structure Refinement of Skutterudite-Type Cobalt Triantimonide, CoSb3. Acta Crystallogr. Sect. C Cryst. Struct. Commun. 1987, 43, 1678– 1679.

(33)

Yamanaka, S.; Enishi, E.; Fukuoka, H.; Yasukawa, M. High-Pressure Synthesis of a New Silicon Clathrate Superconductor, Ba8Si46. Inorg. Chem. 2000, 39, 56–58.

(34)

McWilliams, R. S.; Spaulding, D. K.; Eggert, J. H.; Celliers, P. M.; Hicks, D. G.; Smith, R. F.; Collins, G. W.; Jeanloz, R. Phase Transformations and Metallization of Magnesium Oxide at High Pressure and Temperature. Science 2012, 338, 1330–1333.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(35)

Stackhouse, S.; Stixrude, L.; Karki, B. B. Thermal Conductivity of Periclase (MgO) from First Principles. Phys. Rev. Lett. 2010, 104, 1–4.

(36)

Charvat, F. R.; Kingery, W. D. Thermal Conductivity: XIII, Effect of Microstructure on Conductivity of Single‐Phase Ceramics. J. Am. Ceram. Soc. 1957, 40, 306–315.

(37)

Slack, G. A. Thermal Conductivity of MgO, Al2O3, MgAl2O4, and Fe3O4 Crystals from 3° to 300° K. Phys. Rev. 1962, 126, 427.

(38)

Touloukian, Y. S.; Powell, R.; Ho, C.; Klemens, P. Thermal Conductivity: Nonmetallic Solids, IFI. Plenum 1970.

(39)

Andersson, S.; Bäckström, G. Techniques for Determining Thermal Conductivity and Heat Capacity under Hydrostatic Pressure. Rev. Sci. Instrum. 1986, 57, 1633–1639.

(40)

Dalton, D. A.; Hsieh, W. P.; Hohensee, G. T.; Cahill, D. G.; Goncharov, A. F. Effect of Mass Disorder on the Lattice Thermal Conductivity of MgO Periclase under Pressure. Sci. Rep. 2013, 3, 1–5.

(41)

Chen, B.; Hsieh, W.-P.; Cahill, D. G.; Trinkle, D. R.; Li, J. Thermal Conductivity of Compressed H 2 O to 22 GPa: A Test of the Leibfried-Schlömann Equation. Phys. Rev. B 2011, 83, 132301.

(42)

Imada, S.; Ohta, K.; Yagi, T.; Hirose, K.; Yoshida, H.; Nagahara, H. Measurements of Lattice Thermal Conductivity of MgO to Core‐mantle Boundary Pressures. Geophys. Res. Lett. 2014, 41, 4542–4547.

(43)

Speziale, S.; Zha, C.-S.; Duffy, T. S.; Hemley, R. J.; Mao, H. Quasi-Hydrostatic Compression of Magnesium Oxide to 52 GPa: Implications for the Pressure-VolumeTemperature Equation of State. J. Geophys. Res. Solid Earth 2001, 106, 515–528.

(44)

Chopelas, A. Thermal Expansion, Heat Capacity, and Entropy of MgO at Mantle Pressures. Phys. Chem. Miner. 1990, 17, 142–148.

(45)

Togo, A.; Tanaka, I. First Principles Phonon Calculations in Materials Science. Scr. Mater. 2015, 108, 1–5.

ACS Paragon Plus Environment

Page 20 of 30

Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(46)

Slack, G. A.; Rowe, D. M. CRC Handbook of Thermoelectrics. CRC, Boca Raton, FL 1995, 407.

(47)

Tse, J. S.; Uehara, K.; Rousseau, R.; Ker, A.; Ratcliffe, C. I.; White, M. A.; MacKay, G. Structural Principles and Amorphouslike Thermal Conductivity of Na-Doped Si Clathrates. Phys. Rev. Lett. 2000, 85, 114.

(48)

Tse, J. S.; White, M. A. Origin of Glassy Crystalline Behavior in the Thermal Properties of Clathrate Hydrates: A Thermal Conductivity Study of Tetrahydrofuran Hydrate. J. Phys. Chem. 1988, 92, 5006–5011.

(49)

Pailhès, S.; Euchner, H.; Giordano, V. M.; Debord, R.; Assy, A.; Gomès, S.; Bosak, A.; Machon, D.; Paschen, S.; de Boissieu, M. Localization of Propagative Phonons in a Perfectly Crystalline Solid. Phys. Rev. Lett. 2014, 113, 25506.

(50)

Schopf, D.; Euchner, H.; Trebin, H.-R. Effective Potentials for Simulations of the Thermal Conductivity of Type-I Semiconductor Clathrate Systems. Phys. Rev. B 2014, 89, 214306.

(51)

Tse, J. S.; Ratcliffe, C. I.; Powell, B. M.; Sears, V. F.; Handa, Y. P. Rotational and Translational Motions of Trapped Methane. Incoherent Inelastic Neutron Scattering of Methane Hydrate. J. Phys. Chem. A 1997, 101, 4491–4495.

(52)

English, N. J.; John, S. T. Mechanisms for Thermal Conduction in Methane Hydrate. Phys. Rev. Lett. 2009, 103, 15901.

(53)

Rull-Bravo, M.; Moure, A.; Fernandez, J. F.; Martin-Gonzalez, M. RSC Adv. 2015, 5, 41653.

(54)

Caillat, T.; Borshchevsky, A.; Fleurial, J. Properties of Single Crystalline Semiconducting CoSb3. J. Appl. Phys. 1996, 80, 4442–4449.

(55)

Morelli, D. T.; Caillat, T.; Fleurial, J.-P.; Borshchevsky, A.; Vandersande, J.; Chen, B.; Uher, C. Low-Temperature Transport Properties of P-Type CoSb 3. Phys. Rev. B 1995, 51, 9622.

(56)

Guo, R.; Wang, X.; Huang, B. Thermal Conductivity of Skutterudite CoSb3 from First Principles: Substitution and Nanoengineering Effects. Sci. Rep. 2015, 5,7806 .

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(57)

Zhang, Q.; Li, X.; Kang, Y.; Zhang, L.; Yu, D.; He, J.; Liu, Z.; Tian, Y.;Xu, B. High Pressure Synthesis of Te-Doped CoSb3 with Enhanced Thermoelectric Performance. J. Mater. Sci. Mater. Electron. 2015, 26, 385–391.

(58) Che, J; Cagrin, T.; Deng, W.; Goddard III, W.A., Thermal conductivity of diamond and related materials from molecular dynamics simulations., J. Chem. Phys., 2000, 113, 6888. (59) Turneym J.E.; McGaughey, A.J.H.; Amon, C.H., Assessing the applicability of qunatum corrections to classical thermal conductivity predictions., Phys. Rev. B, 2009, 79, 224305. (60) Baldwin, T. O., Tompson, C. W., X‐Ray Characteristic Temperatures of Some II—VI Ionic Compounds., J. Chem. Phys., 1964, 41, 1420. (61) Tanigaki, K., Shimizu, T., Itoh, K.M., Teraoka, J., Moritomo Y., Yamanaka, S., Mechanism of superconductivity in the polyhedral-network compound Ba8Si46., Nature Mat., 2003, 2, 653. (62) Caillat, T., Borshchevsky, Fleuria, J-P., Preparation and Thermoelectric Properties of pand n-type CoSb3., AIP Conference Proceedings, 1994, 316, 58. (63)

Skylaris, C.-K.; Haynes, P. D.; Mostofi, A. A.; Payne, M. C. Introducing ONETEP: Linear-Scaling Density Functional Simulations on Parallel Computers. J. Chem. Phys. 2005, 122, 84119.

ACS Paragon Plus Environment

Page 22 of 30

Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fig. 1. Comparison of theoretical and experimental thermal conductivity of MgO as a function of temperature at ambient pressure.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 30

Fig. 2. Comparison of theoretical and experimental thermal conductivity of MgO as a function of pressure and 300 K. The experimental data are represented by open triangles with the filled red triangles corrected for the heat capacity (see text).

The dot-dash line is

computed from the theoretical calculated thermal conductivities using cubic anharmonic force constants fitted to an analytical function (cf. eqn. 1 of ref. 8). The MD results were obtained in this study.

ACS Paragon Plus Environment

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fig. 3. Comparison of the heat capacity at constant pressure (Cp) of MgO from theoretical quasiharmonic approximation (QHA) (red line) and estimated from the experimental equation of states (black line) and a second-order polynomial fit (blue line) reported in the Supplementary Materials of ref. 35.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 4. Thermal conductivity of types I (Si46) and II (Si132) empty Si clathrates15 and lattice thermal conductivity of Ba8Si46 (this work); the dotted blue line is aid to the eye. The red star denotes the lattice thermal conductivity of Ba8Si46 calculated from Green-Kubo heatflux autocorrelation formula of a molecular-dynamics calculation wherein the interatomic potential has been fitted to first-principles calculations. The experimental data for Ba8Si46 is also shown.

ACS Paragon Plus Environment

Page 26 of 30

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fig. 5. Theoretical and experimental thermal conductivity of pure CoSb3. The solid blue and dotted lines are the results of harmonic and anharmonic lattice-dynamics calculations using the LDA and GGA density-functional approaches, respectively. The red open circles are the present MD results. The lattice thermal conductivity for a Te-doped CoSb3 at 80 and 300 K are shown as open red diamonds. The blue arrow highlights the thermal conductivity plateau at high temperature. The estimated errors for the calculated thermal conductivity of pure CoSb3 are given in text.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 6. Comparison of theoretical calculated lattice thermal conductivity of Co4Sb12-xTex (x=0.2) with experimental data with x = 0.1 and 0.3 from ref. 52. The estimated errors for the calculated thermal conductivity are given in text.

ACS Paragon Plus Environment

Page 28 of 30

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1. Comparison of theoretical and experimental thermal conductivity (W/m/K) of MgO at high temperature (ambient pressure) and high pressure (300 K).

T (K)

This study

Ref. 8

Expt.

100

152.3 ± 20.6

120.5

300

70.3 ± 8.9

66.0

55.0

600

34.6 ± 4.5

30.1

27.5

900

18.6 ± 2.7

19.1

16.7

1200

13.7 ± 2.0

13.4

12.3

0

70.3 ± 8.7

66

55.0

20

113.5 ± 18.2

129

89.8

47

178.6 ± 25.4

196

137.0

65

211.3 ± 29.8

241

83

253.5 ± 36.6

280

137

357.3 ± 41.0

377

P (GPa)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC graphics

ACS Paragon Plus Environment

Page 30 of 30