Lattice thermal conductivities of two SiO$_2$ polymorphs by first-principles calculation and phonon Boltzmann transport equation

Lattice thermal conductivities of two SiO$_2$ polymorphs, i.e., $\alpha$-quartz (low) and $\alpha$-cristobalite (low), were studied using first-principles anharmonic phonon calculation and linearized phonon Boltzmann transport equation. Although $\alpha$-quartz and $\alpha$-cristobalite have similar phonon densities of states, phonon frequency dependencies of phonon group velocities and lifetimes are dissimilar, which results in largely different anisotropies of the lattice thermal conductivities. For $\alpha$-quartz and $\alpha$-cristobalite, distributions of the phonon lifetimes effective to determine the lattice thermal conductivities are well described by energy and momentum conservations of three phonon scatterings weighted by phonon occupation numbers and one parameter that represents the phonon-phonon interaction strengths.

SiO 2 exhibits many polymorphs including α-quartz and α-cristobalite whose crystal structures are shown in Fig. 1. The numbers of atoms in the unit cells (n a ) are 9 and 12, respectively. Both are made of SiO 4 tetrahedra connected by their vertices. Si atom is located at the center of each tetrahedron and O atoms are at the vertices. The difference of these crystal structures is described by the patterns of the tetrahedron linkages. SiO 4 tetrahedra are more densely packed in α-quartz. As a result, the volume per formula unit is more than ten percent smaller in α-quartz. Their lattice parameters 20,21 are shown in Table. I. Their space-group types are P 3 2 21 (trigonal) for α-quartz and P 4 1 2 1 2 (tetragonal) for α-cristobalite. Both of them in principle have anisotropic thermal conductivity tensors with two independent elements, κ xx and κ zz . Although κ xx and κ zz of α-quartz were reported, 22 only its average value is known for α-cristobalite. 23 The aim of this study is to understand their difference in lattice thermal conductivity. Indeed α-quartz shows much larger anisotropy in lattice thermal conductivity than α-cristobalite as presented in this study. This was investigated from microscopic properties due to phonons. By the long range interaction among atoms and softer low frequency phonon modes, we were required to conduct more careful calculations than that usually we do in conjunction with our software development. 24 thermal conductivity calculations and their analysis are presented in Sec. III. We show similarity and dissimilarity between α-quartz and α-cristobalite in densities of lattice thermal conductivities and distributions of phonon properties as a function of phonon frequency. Then the characteristics of three phonon scatterings are discussed.

A. Computational details
We solved linearized phonon Boltzmann transport equation with single-mode relaxation time approximation. 15,26 We abbreviate this approximation as RTA. For the phonon and lattice thermal conductivity calculations, we employed phonopy 27 and arXiv:1804.01729v2 [cond-mat.mtrl-sci] 9 Apr 2018 phono3py 15 software packages. Unless specially denoted, q-point sampling meshes of 19 × 19 × 19 and 19 × 19 × 14 were used for the lattice thermal conductivity calculations of α-quartz and α-cristobalite, respectively. The isotope scattering effect calculated by the second-order perturbation theory 15,28 was found negligibly small. Therefore it was not included.
The experimental lattice parameters of α-quartz 20 and α-cristobalite 21 were used for all calculations in Sec. III. Choice of the lattice parameters can have a large impact to the lattice thermal conductivity since it is known that decreasing (increasing) lattice parameters increases (decreases) lattice thermal conductivity as has been well studied as pressure dependence of lattice thermal conductivity for many years. [29][30][31][32][33] In Sec. II B, we present calculated lattice thermal conductivity values obtained using experimental and calculated lattice parameters.
Second-and third-order force constants were calculated using the supercell approach with finite atomic displacements of 0.03Å. 15,34 The supercells of 6×6×6 (1944 atoms) and 4×4×4 (768 atoms) of the unit cells were used for the calculations of the second-order force constants of α-quartz and α-cristobalite, respectively. Use of larger supercells is in general important to compute phononphonon scattering channels with better accuracy. For α-quartz, it was necessary to take into account the longrange interaction to remove imaginary acoustic modes in the vicinity of Γ-point. We expect real-space interaction range among three atoms effective for lattice thermal conductivity is relatively shorter than that of the second-order force constants. Therefore, for the thirdorder force constants, we chose 2 × 2 × 2 supercells (72 and 96 atoms). Our supercell choices for α-quartz and α-cristobalite are considered reasonable after the examinations as presented in Sec. II B.
Running many supercell first-principles calculations for the third-order force constants is the most computationally demanding part throughout the lattice thermal conductivity calculation. To omit the computations of parts of force constants in some means, e.g., using realspace cutoff distance, can ease its total computational demand. However we filled all elements of the supercell force constants. Nevertheless our attempts and remarks on using the cutoff distance for computing third-order force constants, that we avoided, are presented in Appendix.
Non-analytical term correction [35][36][37] was applied to dynamical matrices to treat long range dipole-dipole interactions. Though impact of non-analytical term correction to lattice thermal conductivity is often negligible for crystals containing number of atoms in their unit cells such as α-quartz (9 atoms) and α-cristobalite (12 atoms), it turned out to be useful for α-quartz to remove imaginary acoustic modes near Γ-point in conjunction with using the larger supercell.
For the first-principles calculations, we employed the plane-wave basis projector augmented wave method 38 within the framework of density functional theory (DFT) as implemented in the VASP code. [39][40][41] The generalized gradient approximation of Perdew, Burke, and Ernzerhof revised for solids (PBEsol) 42 was used as the exchange correlation potential. A plane-wave energy cutoff of 520 eV was employed. The radial cutoffs of the PAW datasets of Si and O were 1.90 and 1.52Å, respectively. The 3s and 3p electrons for Si and the 2s and 2p electrons for O were treated as valence and the remaining electrons were kept frozen. Reciprocal spaces of the α-quartz supercells used for the calculations of the third-and secondorder force constants were sampled by the 3 × 3 × 3 mesh and at only Γ-point, respectively. The former mesh was shifted by a half grid distance in c * direction from the Γ-point centered mesh. For the α-cristobalite supercells, the reciprocal spaces were sampled by the 2 × 2 × 2 and 1×1×1 meshes with half grid shifts along all three directions from the Γ-point centered meshes, respectively. To obtain atomic forces, the total energies were minimized until the energy convergences became less than 10 −8 eV.
Static dielectric constant tensors and Born effective charge tensors were calculated from density functional perturbation theory as implemented in the VASP code. 43,44 These tensors were symmetrized by their space-group and crystallographic-point-group operations. A sum rule was applied to the Born effective charge tensors following Ref. 37. For these calculations, the plane-wave cutoff energy of 600 eV was used. The reciprocal spaces of the α-quartz and α-cristobalite were sampled by the 12 × 12 × 12 and 8 × 8 × 8 k-point sampling meshes, respectively. The former mesh was shifted by a half grid distance along c * direction and the later mesh was shifted by half grid distances along all three directions from the Γ-point centered meshes.

B. Choices of exchange correlation potentials and convergence criteria
We performed series of lattice thermal conductivity calculations against different exchange correlation potentials, solutions of linearized phonon Boltzmann transport equation, and convergence criteria. We present our calculation results on them. After these examinations, we chose the calculation settings described in Sec. II A, which are considered to give results accurate enough for our discussion.
In Table I, the experimental lattice parameters 20,21 of α-quartz and α-cristobalite and those optimized by calculations are presented. For the calculations, we employed the exchange correlation potentials of PBEsol and local density approximation (LDA). 45 Thermal expansion was not considered in the calculations. Although the calculations show good agreements with the experimental values, we can find that those with PBEsol and LDA slightly overestimate and underestimate the experimental values, respectively.
In Fig. 2, convergences of lattice thermal conductivities with respect to the number of sampling phonon modes  in Brillouin zones are presented. For both of α-quartz and α-cristobalite, the lattice thermal conductivities converge well by ∼10 5 phonon-mode sampling points. Since we needed more sampling phonon modes to converge the curve shapes of spectrum-like plots such as phonon density of states (DOS), we chose the 19 × 19 × 19 and 19 × 19 × 14 q-point sampling meshes for α-quartz and α-cristobalite, respectively. In Table II, experimental and calculated lattice thermal conductivities are presented. For the calculations, we employed RTA 15,26 and direct 11,15 solutions of linearized phonon Boltzmann transport equation, and the obtained values were close each other for α-quartz and α-cristobalite. Therefore, we decided to use the RTA solution, since, compared with the direct solution, it has an advantage in analyzing results more easily and intuitively by its closed form of lattice thermal conductivity formula. Due to crystal symmetries of α-quartz (trigonal) and α-cristobalite (tetragonal), both lattice thermal conductivity tensors have only two degrees of freedom, κ xx and κ zz . α-quartz exhibits largely anisotropic lattice thermal conductivity whereas that of α-cristobalite is more isotropic as shown in Table II. From the experimental measurement of α-quartz by Kanamori et al., 22 the ratio κ zz /κ xx is around 2, which is well reproduced by our calculation. However each of κ xx and κ zz from the calculation underestimates the experimental values. There is another experimental measurement of powder sample reported by Kunugi et al. 23 By taking κ av = (2κ xx +κ zz )/3 as an averaged value along orientations, the calculated value is found to be close to the experiment. In the same report by Kunugi et al., they also showed the measurement of powder α-cristobalite, which agrees well with the averaged value by the present calculation. TABLE III. Calculated lattice thermal conductivities κ (W/m-K) of α-quartz and α-cristobalite at 300 K with respect to the choices of lattice parameters (see Table I It is not always the case that we can fortunately refer to experimental lattice parameters on lattice thermal conductivity calculations. Therefore it is of interest to see how much different lattice thermal conductivities are calculated using the lattice parameters determined by the first-principles calculations and those calculated with the experimental lattice parameters. The results are shown in Table III. For the calculated lattice parameters by PBEsol (LDA) overestimated (underestimated) the experimental lattice parameters as shown in Table I, smaller (larger) lattice thermal conductivities were obtained following the general trend of the volume dependence. When using the same experimental lattice parameters, the lattice thermal conductivities calculated with PBEsol were obtained larger than those with LDA for both α-quartz and α-cristobalite. From these calculations, we can see distinguishable effects by the choices of the exchange correlation potentials: one is in determining lattice parameters and the other is in calculating forces on atoms. However since the values and the ratios κ zz /κ xx in Table III are close enough, any choice given here is found a reasonable choice unless we expect too good quantitative agreements between calculations and experiments. We investigated the effects on calculated lattice thermal conductivities by different choices of supercell size used for the calculation of the third-order force constants and finite atomic displacement distance and plane-wave cutoff energy used for the calculations of the second-and third-order force constants. The k-points of the supercell reciprocal spaces were sampled with equivalent density meshes to those written in Sec. II A except for that of 3 × 3 × 2 supercell of α-cristobalite where the 2 × 2 × 2 sampling mesh shifted in half grid distances along all directions from the Γ-point centered mesh was used. These results show, for both α-quartz and α-cristobalite, that 2 × 2 × 2 supercells are the reasonable choices considering the tradeoff of convergences of the lattice thermal conductivity values and the required computational demands (see Appendix) with respect to our current computational resource. It also shows the use of the unit cells for third-order force constants calculations is not a bad choice if a purpose is the rough estimation.
The choice of 0.05Å displacement distance induces decrease of lattice thermal conductivity for α-quartz. This is considered due to inclusion of higher order anharmonicity. In general, decreasing the displacement distance, numerical error in force constants is magnified. The results by the choice of 0.01Å displacement distance give similar results with those by 0.03Å. This means that the numerical errors and inclusions of higher order anharmonicity are managed to be small by the choice of 0.03Å displacement distance for our computer simulation settings.

III. RESULTS AND DISCUSSION
In RTA, lattice thermal conductivity κ is written in a closed form: 26 where N and V 0 are the number of unit cells in the system and volume of the unit cell, respectively. The suffix λ represents the phonon mode as the pair of phonon wave vector q and branch j, λ ≡ (q, j), and similarly we denote −λ ≡ (−q, j). C λ is the mode heat capacity given as where ω λ = ω(q, j) is the phonon frequency. v λ is the phonon group velocity defined as the gradient of the phonon energy surface: v λ = ∇ q ω(q, j).
τ λ is the single-mode relaxation time and we use phonon lifetime as τ λ . We calculated phonon lifetime where with n λ = [exp( ω λ /k B T ) − 1] −1 as the phonon occupation number at equilibrium. In these equations, and k B denote the reduced Planck constant and Boltzmann constant, respectively. Φ λλ λ gives the phononphonon interaction strength among three phonons calculated from second-and third-order force constants. ∆(q + q + q ) ≡ 1 if q + q + q = G otherwise 0, where G is the reciprocal lattice vector. This constraint comes from the lattice translational invariance that appears inside Φ λλ λ , 15 however it is made visible in Eq. (4) for the analysis given below. More details such as the phase convention, coefficients, and Φ λλ λ are found in Ref. 15. Phonon band structures and DOS of α-quartz and αcristobalite are shown in Figs. 3 (a) and (b), respectively. These phonon structures in their shapes show reasonable agreements with previous calculations and experiments reported in Refs. 36, 47-51. Between α-quartz and α-cristobalite, their total and partial DOS curves are analogous. In detail, the position of the first peak of α-quartz from 0 THz is located at higher phonon frequency than that of α-cristobalite. Their first peak positions roughly correspond to M-and L-points of α-quartz and M-point of α-cristobalite in respective phonon band structures. These low phonon modes are considered to be made of rigid unit motions of SiO 4 tetrahedra, [52][53][54] i.e., the phonon band structures at low frequencies reflect the different styles of the tetrahedron linkages.
To visualize phonon mode contribution to lattice thermal conductivity, we define densities of lattice thermal conductivities κ(ω) as Compared with phonon DOS written as 1/N λ δ(ω − ω λ ), Eq. (7) is considered as a weighted DOS and each weight C λ v λ ⊗v λ τ λ /V 0 is understood as microscopic contribution of phonon mode λ to lattice thermal conductivity at ω λ . In Fig. 4 (a), κ(ω) of α-quartz and αcristobalite are drawn as a function of phonon frequency. We can see large peaks below 5 THz, where the phonon modes determine more than halves of κ xx and κ zz of α-quartz and α-cristobalite. The curve shapes of κ(ω) are similar to those of the phonon DOS below their first peaks. Therefore it is considered that the number of states is the dominating factor of the lattice thermal conductivities in these phonon frequency ranges. Above 5 THz, κ(ω) are relatively small, however they contribute little by little to κ up to ∼25 THz. Anisotropy of lattice thermal conductivity, i.e., the ratio κ zz /κ xx , is larger in α-quartz than in α-cristobalite. The phonon mode contributions to the anisotropic κ are discussed using cumulative lattice thermal conductivity given by Obviously lim ω→∞ κ c (ω) = κ from Eq. (6). The ratios κ c zz (ω)/κ c xx (ω) are shown in Fig. 5, where α-quartz and α-cristobalite present similar behaviours, although their intensities are different. Increasing phonon frequency from 0 THz, their ratios rapidly increase at low phonon frequencies and start to decrease gently until the ratios become κ zz /κ xx . This difference is exhibited in distributions of v λ ⊗ v λ that are written in analogy to κ(ω) of Eq. (7) as w(ω) are shown in Fig. 4 (b). Below 5 THz, the ratio between w zz (ω) and w xx (ω) is clearly larger in α-quartz than in α-cristobalite. Comparing Figs. 4 (a) and (b), increasing phonon frequency, κ(ω) more quickly decrease after first large peaks than w(ω) in both α-quartz and α-cristobalite. This is due to phonon frequency dependencies of C λ τ λ , however the effect of C λ to the curve shapes of κ(ω) with respect to those of w(ω) is small since C λ is approximately constant ∼k B at 300 K below 10 Hz. In Fig. 6, τ λ are plotted by dots as a function of phonon frequency. At lower phonon frequencies, phonons tend to have longer lifetimes and decrease quickly their lifetimes increasing phonon frequency below 5 THz. Both of α-quartz and αcristobalite show the same trend but with different rate of decrease, which clearly impacts to the shapes of κ(ω) in Fig. 4 (a), e.g., κ(ω) of α-cristobalite corresponding to the second peak of w(ω) at ∼3 THz is removed by the decrease of τ λ . Recalling Eq. (4), τ λ is constructed from the wave vector constraint ∆(q + q + q ), weighted energy conservation N λ λ (ω λ ), and |Φ λλ λ | 2 . To make our discussion simple, we replace |Φ λλ λ | 2 by a constant valueP if q + q + q = G or by 0 otherwise. As an attempt, we useP =P av defined as an average of |Φ λλ λ | 2 bỹ where P λ is that for one phonon mode: 15 Since (3n a ) 2 NP av of α-quartz and α-cristobalite give the equivalent values as shown in Table.V, we expect that they have similar phonon-phonon interaction strengths.
WithP , Γ λ (ω) is reduced tõ In Eq. (12), the summation on the right hand side is made of three phonon scattering channels weighted by phonon occupation numbers, which can be computed from the second-order force constants. The lattice thermal conductivities calculated withP =P av andτ λ = (2Γ λ ) −1 , that we denoteκ, are presented in Table V. These values are one order of magnitude smaller than the values in Table II, however the anisotropiesκ zz /κ xx are well reproduced, and as shown in Fig. 7, the curve shapes of the densities of lattice thermal conductivities, denoted byκ(ω), are almost identical to those of κ(ω) presented in Fig. 4 (a). In Fig. 8, P λ of α-quartz and α-cristobalite are plotted as a function of phonon frequency. Their distributions are similar except at low phonon frequency do- . Between 0 to 15 THz, P λ are roughly constant, by which, apart from their different magnitudes, similar phonon frequency dependencies ofτ λ to those of τ λ are obtained as shown in Fig. 9. This enables the curve shapes ofκ(ω) to become equivalent to those of κ(ω). Since more than 90% of the lattice thermal conductivities of α-quartz and α-cristobalite are recovered in κ(ω) below 15 THz, having a good estimate of the constant value, e.g.,P ∼P av × 10 −1 , it is considered possible to predict the lattice thermal conductivities without computing third order force constants. P λ start to increase from ∼15 THz to the phonon band gap at ∼25 THz. The two small domains above 30 THz correspond to respective two localized phonon bands. The ratio of Si and O partial DOS gradually increases by increasing phonon frequency below 15 THz. This represents that SiO 4 rigid units vibrate translationally at lower phonon frequencies and rotationally at increasing .τ λ , phonon lifetimes of α-quartz and α-cristobalite calculated withP =Pav at 300 K as a function of phonon frequency (black dots). To compare, τ λ × 10 −1 (Fig. 6) are shown as the grey dots behind the black dots. The points are sampled on the 19×19×19 mesh for α-quartz and 19×19×14 mesh for α-cristobalite in the respective Brillouin zones.
phonon frequencies. Above 15 THz, it is considered that the larger P λ , i.e., larger anharmonicity, arises due to phonons that distort SiO 4 tetrahedron units.

IV. SUMMARY
The lattice thermal conductivity calculations were performed for α-quartz and α-cristobalite using firstprinciples anharmonic phonon calculation and linearized phonon Boltzmann transport equation. Since direct and RTA solutions gave similar values of the lattice thermal conductivities that also agree well with the experimental values, we focused on our discussion using the RTA solutions and phonon frequency dependencies of the phonon properties. The densities of the lattice thermal conductivities κ(ω) show the characteristic differences of phonon mode contributions to the lattice thermal conductivities between α-quartz and α-cristobalite. Below 2 THz for α-cristobalite and 3 THz for α-quartz, phonon DOS and v λ ⊗ v λ determines the shapes of κ(ω). Above 5 THz, κ(ω) becomes much smaller than those below 5 THz following the phonon frequency dependence of τ λ . The large difference of anisotropies in the lattice thermal conductivities of α-quartz and α-cristobalite was found. This is mainly attributed by the distributions of the phonon group velocities below 5 THz. The distributions of the phonon lifetimes effective to determine the lattice thermal conductivities around room temperature were well described by the momentum conservation ∆(q+q +q ), the energy conservation weighted by the phonon occupation numbers N λ λ (ω λ ), and a parameterP that represents the phonon-phonon interaction strengths. Appendix A: Effect of using real-space cutoff to calculate supercell third-order force constants Use of real-space cutoff distance to compute thirdorder force constants in the supercell approach may drastically reduce computational demand of lattice thermal conductivity calculation. However it should be used carefully since the side effect such as degradation of the numerical quality has not been well understood. In this Appendix, we provide our examinations on the effect of using a cutoff distance for the third-order force constants calculations. There are many possible ways to cut off third-order force constants. Below, we explain our scheme and show the convergence analysis.

Scheme
We calculate supercell third-order force constant element from two atomic displacements and a force on an atom by, 15 where u α (lκ) means the finite displacement of the atom at the position r(lκ) along α-th Cartesian axis. The indices l and κ denote the lattice point and the atom in the unit cell, respectively. F γ [l κ ; u(lκ), u(l κ )] gives the force that the atom l κ experiences by two atomic displacements u(lκ) and u(l κ ). Here it is assumed that we can obtain forces on all atoms in the supercell at once by each supercell calculation with a pair of atomic displacements. This assumption is currently normal in the DFT calculations since the computation of forces from existing electronic wave function requires relatively small computation. Our cutoff distance R cut is used to collect all the displaced atomic pairs whose distances |r(lκ) − r(l κ )| 2 are shorter than R cut . The set of these pair displacements fills all supercell third-order force constant elements except for the elements whose three atoms are mutually more distant than R cut .

Results
In this section, we present calculated lattice thermal conductivities using different cutoff distances and see the convergences in α-quartz and α-cristobalite using different supercell sizes. We employed 3×3×2 and 2×2×2 supercells and unit cells for these examinations. The computations of third-order force constants using the 3×3×2 supercells were computationally very demanding for our current computational resource to fill all the elements, but not with the 2 × 2 × 2. For α-quartz, one 3 × 3 × 2 supercell calculation was five times more computationally demanding than one 2 × 2 × 2 supercell calculation. For α-cristobalite, that was nine times because of the denser k-point sampling for the 3 × 3 × 2 supercell calculation of α-cristobalite.
The purpose to use the cutoff distance is to obtain accurate third-order force constants with reasonable computational demand though it is safer to compute all elements of supercell force constants to avoid sudden cut of those elements since it is difficult to predict what happens after Fourier transformation of the third-order force constants with the cut.
In Figs. 10 (a), (b), and (c), the lattice thermal conductivities of α-quartz calculated against the cutoff distances are shown for three different supercell sizes. The lattice thermal conductivities generally decrease increasing the cutoff distance in these supercell sizes. It looks that each lattice thermal conductivity converges toward its rightmost point that corresponds to the full calculation where all elements of the supercell third-order force constants were filled. In Fig. 10 (b), at the rightmost point, the lattice thermal conductivity increases in contradiction to the tendency of decreasing with increasing the cutoff distance. This is considered a visible effect of the cut of the supercell third-order force constants elements. For αcristobalite as shown in Figs. 10 (d), (e), and (f), the convergence is achieved at relatively shorter cutoff distance of ∼4Å. This is about the distance between two atoms in neighboring SiO 4 tetrahedra. However the calculation of the third-order force constants with the 3 × 3 × 2 supercell using ∼4Å cutoff distance is already more computationally demanding than the full calculation with the 2 × 2 × 2 supercell. Therefore the supercell size has to be chosen systematically along with the choice of the cutoff Lattice thermal conductivities at 300 K with respect to cutoff distances of atomic pairs used to compute third-order force constants employing (a) α-quartz 3 × 3 × 2 supercell (162 atoms), (b) α-quartz 2 × 2 × 2 supercell (72 atoms), (c) α-quartz unit cell (9 atoms), (d) α-cristobalite 3 × 3 × 2 supercell (216 atoms), (e) α-cristobalite 2 × 2 × 2 supercell (96 atoms), and (f) α-cristobalite unit cell (12 atoms). The selected cutoff distances are those closest to but below (a) 2, . . . , 11Å, (b) 2, . . . , 8Å, and (c) 2, 3, 3.5, and 4Å, (d) 2, . . . , 12Å, (e) 2, . . . , 10Å, and (f) 2, . . . , 5Å, respectively. The filled circles depict κzz/2 for α-quartz and κzz for α-cristobalite, and the open circles show κxx. The cross symbols present the numbers of supercells with displacements that were used to compute the third-order force constants with the respective cutoff distances. The rightmost points correspond to the results obtained without using the cutoff distances. Lines are eye guides.
distance. Comparing Figs. 10 (a) and (d), we can see lattice thermal conductivity of α-cristobalite converges more quickly than that of α-quartz. For α-quartz, it is difficult to define the convergence criterion to choose the cutoff distance for the accurate lattice thermal conductivity calculation.
For a purpose of the rough estimation, any choice of the cutoff distance and supercell size seems acceptable. In the case of α-quartz and α-cristobalite, the first nearest neighbor distance is well isolated because of SiO 4 tetrahedra. This may be the reason why the atomic interaction range effective to determine lattice thermal conductivity is found to be short.