Nanoscale thermal transport at PbTe grain boundaries: a neural-network-potential molecular dynamics study

Abstract

Grain boundaries (GBs) are nano- and mesoscale defects that play a key role in tailoring lattice thermal conductivity κ_l and enhancing thermoelectric conversion efficiency. However, even in the prototypical thermoelectric material PbTe, the relationship between GB atomic structure and lattice thermal transport remains poorly understood. Here, we combine a neural-network potential (NNP) with perturbed molecular dynamics to determine κ_l in the directions normal and parallel to the GB planes for 24 symmetric tilt GBs in PbTe. The NNP accurately reproduces GB energetics and finite-temperature atomic forces, enabling large-scale simulations of GB models containing thousands of atoms. Along the GB-normal direction, κ_l shows a tendency to decrease with coordination-number change, and the corresponding GB thermal resistances vary by a factor of seven depending on the GB structure. Along the GB-parallel directions, κ_l is influenced by bond-length distortion and the crystallographic orientation of grains, and the effect of GBs extends over a spatial range of about 10 nm from the GB planes. These results suggest that thermal transport parallel to GB planes is also non-negligible, particularly for nanocrystalline samples. Atom-resolved thermal conductivities exhibit only weak correlations with local structural descriptors, highlighting the necessity of macroscopic descriptors that encode the GB network and its relation to the direction of heat flux. These findings provide fundamental insights into GB-dependent phonon transport in PbTe and offer a basis for data-driven GB engineering and microstructure design to further improve thermoelectric performance.

---

1 Introduction

Thermoelectric devices can directly convert waste heat into electricity and are therefore regarded as an environmentally friendly technology.^1,2 The performance of a thermoelectric material is characterized by the dimensionless figure of merit ZT = S²σT/(κ_l + κ_e), where T is the absolute temperature, S is the Seebeck coefficient, σ is the electrical conductivity, and κ_l and κ_e are the lattice and electronic thermal conductivities. Because S, σ, and κ_e are all functions of the carrier concentration and are strongly interdependent, an important strategy for improving thermoelectric conversion efficiency is to reduce κ_l while maintaining favorable electronic transport properties.³

Lead telluride (PbTe) is a representative thermoelectric material that exhibits high conversion efficiency in the intermediate temperature range.^4,5 The highest ZT of 2.8 at 850 K has been reported for p-type PbTe,⁶ and in recent years high ZT values of ∼2 have also been achieved in n-type PbTe.⁷ Despite the highly symmetric rocksalt structure, pristine PbTe shows a low thermal conductivity of approximately 2.3 W m⁻¹ K⁻¹ at 300 K,^8–10 due to its heavy constituent elements and strong anharmonicity.^11–13 To further reduce κ_l, various nanostructuring strategies have been extensively employed for PbTe.^3,14,15 These include point defects,^6,16 nanoprecipitates,^17,18 grain boundaries (GBs),^19,20 and dislocations,^21,22 which together create hierarchical defects over multiple length scales for phonon scattering. This strategy is effective for PbTe because its electronic transport properties are relatively insensitive to defect scattering compared with those of other thermoelectric materials such as Mg₃Sb₂.²³

GBs are ubiquitous in polycrystalline materials and are commonly present in PbTe as nano- and mesoscale phonon-scattering sources.^3,4 Reducing the grain sizes down to 20–100 nm can lower κ_l to 1.7 W m⁻¹ K⁻¹.¹⁹ In many previous studies, the effect of GBs has been treated in an averaged manner, and their influence has been discussed primarily in terms of GB populations (grain size). However, GBs possess a vast variety of possible structures and associated properties, which originate from large configurational freedom depending on grain misorientation.^24–26 For example, Isotta et al. measured the local thermal conductivity in polycrystalline rocksalt SnTe, and demonstrated that the thermal resistances of individual GBs vary depending on the misorientation angle or GB atomic structure by a factor of 2.5.²⁴ Wu et al. measured the mobility of electronic carriers across six GBs in PbTe and showed that the potential barrier height increases at high-angle GBs.²⁵ These findings indicate that, for the design and control of thermoelectric performance via GBs, it is crucial to reveal the relationship between the GB atomic structure and the corresponding transport properties.

Atomistic simulations, such as density functional theory (DFT) calculations, provide a powerful approach for treating the large configurational freedom of GBs and for clarifying correlations between GB atomic structure and thermal conductivity.²⁷ However, the direct integration of DFT calculations, GB modeling, and thermal conductivity calculations is computationally demanding. Alternatively, many previous studies have used empirical interatomic potentials with lower computational cost to assess qualitative trends in GB structures and κ_l, as reported for Si,²⁸ MgO,^29–31 ZnO,³² and SrTiO₃.³³ Recently, machine-learning interatomic potentials (MLIPs) trained on DFT data have been developed to retain DFT accuracy while being several orders of magnitude faster than DFT calculations.³⁴ This development has made it possible to perform large-scale thermal conductivity calculations, such as molecular dynamics (MD) simulations with thousands of atoms and several million time steps.^35,36 For PbTe, Qin et al. evaluated κ_l of twin boundaries in PbTe using an MLIP, as well as its synergistic effect with point defects.³⁷ Nevertheless, applications of MLIPs to GB thermal transport remain limited to few systems, and comprehensive analyses across crystallographically diverse GBs are necessary for establishing structure–property relationships.

Here, we systematically investigated κ_l of PbTe GBs using a neural-network potential (NNP), which was trained on DFT data including several representative GB structures. We modeled 24 symmetric tilt GBs with the [001] and [11̄0] rotation axes and demonstrated that their energetics and atomic forces are accurately predicted by the NNP. For the predicted lowest-energy GB structures, NNP-based MD simulations were performed to determine κ_l and atom-resolved thermal conductivities along the directions normal and parallel to the GB planes. To clarify differences in thermal transport between the two rotation axes, mode-decomposed analyses of the longitudinal and transverse phonon contributions were also carried out. Finally, to identify the structural factors that govern GB thermal conductivity, we examined how κ_l and atom-resolved thermal conductivities correlate with various structural metrics, including coordination-number change and bond-length distortion.

2 Computational methods

2.1 DFT calculation

The Vienna Ab initio Simulation Package (VASP)^38,39 was used to perform DFT calculations within the projector augmented wave (PAW) framework.^40,41 The 5d, 6s and 6p electrons of Pb and the 5s and 5p electrons of Te were treated as valence states in the PAW pseudopotential. The generalized gradient approximation revised for solids (GGA-PBEsol)^42,43 was used for the exchange-correlation functional. The cutoff energy was set to 500 eV, and the convergence criterion for self-consistent electronic-structure calculations was set to 10⁻⁶ eV. A 6 × 6 × 6 grid mesh was used for k-point sampling of the conventional unit cell. Ab initio MD simulations were performed using a Langevin thermostat in the NVT and NPT ensembles, combined with the Parrinello–Rahman dynamics^44,45 in the NPT ensemble.

2.2 Neural network potential

An NNP for PbTe was developed using our in-house code based on the architecture proposed in previous studies,^46–48 with two hidden layers each having 60 nodes. To encode the crystal structures, a structural descriptor constructed from Chebyshev polynomials⁴⁹ was used with 16 radial terms and 48 angular terms for each element. Their cutoff radius was set to 7.5 Å. Further details of the implementation can be found in our recent study.³⁶ An extended Kalman filter^50,51 combined with mini-batch learning was employed to train the NNP.

The training dataset was constructed using not only the perfect crystal but also systems with point defects, surfaces, and GBs as reference structures to cover a wide range of atomic environments (see the SI for details). First, the NNP was trained using atomic structures obtained from DFT calculations in the following procedure. For structural relaxation, initial structures were generated by randomly displacing atoms and distorting the cell dimensions and shape and then relaxed for a few iterations. MD simulations were performed at 100–1600 K for a total duration of 500 fs to generate snapshots every 4 fs. Second, the trained NNP was incorporated into structural relaxation and MD simulations to sample atomic configurations primarily near equilibrium states. The DFT data for the sampled structures were collected by performing single-point DFT calculations and added to the original training dataset. The expanded dataset was used to retrain the NNP, and its predictive accuracy was further improved by repeating the second step.

2.3 Modeling of grain boundaries

The trained NNP was used to identify low-energy atomic structures of the GBs listed in Table 1. As illustrated in Fig. 1(a), the simulation cell of each GB was constructed with three-dimensional periodic boundary conditions, thus containing two crystallographically equivalent GBs. In the Cartesian coordinate system, the x-axis was aligned with the [001] or [11̄0] rotation axis, the y-axis was oriented perpendicular to the rotation axis and parallel to the GB planes, and the z-axis was aligned perpendicular to the GB planes.

Table 1 Symmetric tilt GBs examined for the [001] and [11̄0] systems. Here, θ represents the misorientation angle between two grains

Rotation axis	Σ number	GB plane	θ (°)
[001]	25	(710)	16.26
	13	(510)	22.62
	17	(410)	28.07
	5	(310)	36.87
	29	(520)	43.60
	5	(210)	53.13
	17	(530)	61.93
	13	(320)	67.38
	25	(430)	73.74
	41	(540)	77.32
[11̄0]	51	(551)	16.10
	19	(331)	26.53
	9	(221)	38.94
	11	(332)	50.48
	33	(554)	58.99
	3	(111)	70.53
	17	(334)	86.63
	17	(223)	93.37
	3	(112)	109.47
	33	(225)	121.01
	11	(113)	129.52
	9	(114)	141.06
	27	(115)	148.41
	51	(117)	157.16

---

Fig. 1 Workflow of GB thermal conductivity analysis. (a) Schematic of a symmetric tilt GB. (b) Determination of GB atomic structure by rigid translation. (c) Evaluation of κ_l and κ^atom_l,i and their correlations with GB structure.

For each GB, multiple structures were initially generated by rigidly translating one grain relative to the other in the x and y directions, as illustrated in Fig. 1(b). All structures were then relaxed by allowing the z-axis length and the atomic positions to change, while fixing the x- and y-axis lengths to the corresponding bulk values. The convergence criterion for atomic forces was set to 5 × 10⁻³ eV Å⁻¹. Note that for several GBs with the [11̄0] axis, the lowest-energy structures could not be obtained with structural relaxation alone, as this process explored only local minima with higher GB energies. This often occurs when the true lowest-energy structure significantly deviates from the atomic configuration of the perfect crystal.^52,53 To increase the probability of escaping from local minima, an NNP-based MD simulation was also performed for 500 ps at 400 K for each initial structure generated by the rigid translation described above. In this MD simulation, snapshots were obtained every 10 ps from 310 ps to 500 ps. The lowest-energy structure was identified by applying structural relaxation to all snapshots and was then used to analyze its GB energy, atomic environment and thermal conductivity.

2.4 Perturbed molecular dynamics

κ_l of the lowest-energy GB structures at 300 K were calculated using the PMD method. For these calculations, we employed supercells with a length of approximately 100 Å along the GB-normal direction (z-axis), corresponding to a separation of 50 Å between the two GB planes, and with cell lengths exceeding 45 Å in the GB-parallel directions (Fig. 1(c)).³⁶ In the PMD method,⁵⁴ a perturbation within the linear response regime is applied to the atoms to enhance a heat flux in a specific direction, and κ_l is evaluated from the resulting non-zero time-averaged heat flux. The microscopic heat flux J = (J_x, J_y, J_z) for many-body potentials^55–57 is defined aswhere J_i is the contribution of atom i to J, V is the system volume, m_i and v_i are the mass and velocity of atom i, I is the second-rank unit tensor, U_i is the potential energy of atom i, and r_ij is the position vector from atom i to atom j. When a perturbation is imposed along the x direction, the lattice thermal conductivity in the x direction can be obtained from the following expression:where F_ext is the magnitude of the perturbation, T is the temperature, and t is the time. Since the heat flux can be decomposed into atomic contributions (J_i), the contribution of each atom to κ_l (κ_l,i) can also be evaluated.⁵⁸ However, κ_l,i is an extensive variable that is inversely proportional to the cell volume. Multiplying it by the number of atoms contained in the cell, N, provides a simple way to convert it into an intensive variable, which is referred to as the atomic thermal conductivity.³⁶

κ^atom_l,i = Nκ_l,i. (3)

With this definition, the average of κ^atom_l,i corresponds to κ_l. k^atom_l,i can be mapped to evaluate the spatial dependence of thermal conduction. For example, the local thermal conductivity projected onto the yz plane can be calculated as follows:

where l_x is the cell length in the x direction, which is perpendicular to the projected plane. The Gaussian width parameter σ was set to 0.6 and 2.0 Å.

In addition to the atomic decomposition, the heat flux was decomposed into wave components, which allowed us to extract the longitudinal and transverse contributions to κ_l from the same perturbed MD simulations.³¹ The spectral κ_l (phonon-mode-resolved contributions to κ_l) was also calculated for pristine PbTe by projecting atomic trajectories from perturbed MD simulations onto the phonon modes derived from lattice dynamics calculations.⁵⁸ The detailed procedure for calculating the total and decomposed κ_l using the PMD method is similar to that employed in previous studies^35,36,59,60 and is provided in the SI.

We estimated κ_l normal to the GB plane into the GB thermal resistance, R_GB, using the following expression:⁶¹

R_GB = d(1/κ_GB − 1/κ_perf), (5)

where κ_GB and κ_perf are the lattice thermal conductivities of the GB model and the perfect crystal, respectively, and d is the separation between two GB planes. The calculated κ_perf at 300 K was 2.2 W m⁻¹ K⁻¹, in good agreement with experiments.^8–10

2.5 Grain boundary characterization

To investigate the correlation between κ_l and GB structure, several structural metrics of GBs were evaluated. The GB energy, ΔE_GB, was defined aswhere E_GB and E_perf are the energies of the GB model and the perfect crystal, N_GB and N_perf are the number of atoms contained in the respective simulation cells, and A is the GB cross-sectional area. The division by 2A accounts for the presence of two identical GB planes in the model. Similarly, the excess volume was defined as follows:where V_GB and V_perf are the volumes of the GB model and the perfect crystal. The coordination number of atom i, N_i, was calculated in a Fermi-smeared fashion as follows:^62,63where r_ij is the distance from atom i to its neighbor atom j, r_mid = 0.5(r_1NN + r_2NN), and σ = 0.2(r_2NN − r_1NN), where r_1NN and r_2NN are the first and second nearest neighbor distances in the perfect crystal. Here, only Te atoms are treated as neighbor atoms of Pb, and vice versa. To use N_i as an indicator for the entire GB, we defined the coordination-number change as the sum of the deviations of the coordination numbers from that in the perfect crystal, as follows:where N_perf is the coordination number in the perfect crystal (5.579). The coordination-number change arises primarily from reduced coordination in the vicinity of GBs. Distortion index of atom i, D_i, was defined as follows:where n is the number of neighboring atoms within the cutoff of r_mid. Note that, in the original definition,⁶⁴ the average bond length within the coordination polyhedron is used instead of r_1NN, whereas in this study we use r_1NN in order to quantify deviations from the perfect crystal. We refer to the sum of D_i divided by 2A as the bond-length distortion, which serves as a structural metric of the GB model.

The smooth overlap of atomic positions (SOAP) descriptor was employed to investigate the distributions and differences in local atomic environments in all GB models. The SOAP vectors were calculated using the DScribe code,^65,66 with the cutoffs of r_2NN and 25 Å, 12 radial basis functions, the maximum degree of spherical harmonics of 12, and a Gaussian width of 0.5 Å. Spherical Gaussian-type orbitals were used as the radial basis functions. We applied a distance-dependent weighting of the atomic density using the “poly” option with parameters c = 1 and m = 1.⁶⁷ In addition, the SOAP vectors were compressed in an element-agnostic manner,⁶⁸ since both Pb and Te possess face-centered-cubic sublattices. We defined the difference between the SOAP vectors in the perfect crystal and atom i in the GB model as a local distortion factor, S_i, as follows:³⁰

where p_i and p_perf are the SOAP vectors of atom i in the GB model and an atom in the perfect crystal.

3 Results and discussion

3.1 Grain boundary structure and energy

This section examines the accuracy of the NNP and summarizes the lowest-energy atomic structures. Fig. 2(a) shows that the NNP accurately predicts DFT GB energies, especially for the [001] system, where all data points lie close to the diagonal line with no significant deviation. Importantly, the NNP maintains its predictive power even for GBs that were not included in the training dataset (Table S1). This indicates broad coverage of GB atomic environments without severe overfitting. For the [11̄0] system, the predictive accuracy is lower than that for the [001] system, particularly for high-energy GBs. Nevertheless, the lower-energy structures still lie close to the diagonal line. This suggests that the NNP retains sufficient accuracy to predict the lowest-energy atomic structures without direct DFT calculations for both the [001] and [11̄0] systems.

Fig. 2 Predictive accuracy of the NNP with respect to (a) GB energies and (b) and (c) atomic forces for the GBs with Σ ≤ 17. Results for the [001] and [11̄0] systems are shown on the left and right panels, respectively. For the GB energies, each data point corresponds to the GB energy of an atomic structure that was obtained by applying structural relaxation to an initial configuration generated by the rigid translation, as mentioned in Section 2.3. The DFT value for the relaxed structure was obtained by single-point DFT calculations. Note that each GB has multiple data points associated with different atomic structures, representing not only the lowest-energy structure but also local-minimum structures. To evaluate atomic forces, NNP-based MD simulations were performed in the NVT ensemble at 200, 400 and 600 K. Snapshots generated at each temperature were then used to calculate the NNP and DFT values. In (c), the difference between NNP and DFT forces of each atom is plotted as a function of its coordinate along the z-axis. The gray-shaded area indicates the position of each GB.

The increased errors observed for the [11̄0] system arise from the high-energy structures of the Σ11(332), Σ17(223) and Σ17(334) GBs, which were not included in the training dataset. In future work, it would be desirable to expand the present training dataset by including GBs with various crystallographic characters, with the ultimate goal of developing the NNP generalizable to any GBs.

The NNP accuracy for MD simulations is evaluated by comparing the NNP and DFT values of atomic forces (Fig. 2(b)). Significant deviations are absent at the three temperatures even for the GBs not included in the training dataset. Table S2 shows that the mean absolute errors remain small, ranging 27–60 meV Å⁻¹ across 200–600 K. Fig. 2(c) shows that the average NNP-DFT force difference is smaller than approximately 80 meV Å⁻¹ and 100 meV Å⁻¹ for the [001] and [11̄0] systems, respectively, although it becomes larger near the GBs. Overall, the NNP reproduces atomic forces sufficiently well for quantitative analyses of κ_l and lattice vibrational properties of GBs.

Atoms at the Σ17(334) GB show larger errors than the other GBs for a few snapshots, ranging approximately 120–140 meV Å⁻¹. Nevertheless, its calculated κ_l follows a same trend as the other GBs (see below), supporting acceptable NNP accuracy for this GB.

Fig. 3 displays the atomic structures at 0 K and 300 K for three representative GBs. All atomic structures of the GBs are provided in the SI. As shown in Fig. 3(a), the Σ3(111) GB corresponds to a coherent twin boundary,^69,70 which is accurately reproduced by the NNP. This GB exhibits essentially the same atomic structure at 0 K and 300 K. In contrast, several GBs exhibit distinctly different atomic structures at the two temperatures. For example, the Σ25(430) GB (Fig. 3(b)) has an asymmetric structure and a symmetric structure with respect to the GB plane at 0 K and 300 K, respectively. This is also the case for the Σ51(551) GB (Fig. 3(c)). At 300 K, its atomic structure contains the [11̄0] columns with a lower Pb-atom density (red arrows), as also observed in several [11̄0] GBs.

Fig. 3 Atomic structures of (a) the Σ3(111), (b) the Σ25(430) and (c) the Σ51(551) GBs at 0 K and 300 K. The finite-temperature structures were obtained by equilibrating the simulation cells and subsequently averaging the atomic positions over MD simulations in the NVT ensemble.

Fig. 4 quantifies deviations in atomic structures of the GBs from the perfect crystal by evaluating three metrics, which are correlated with their κ_l in Sections 3.3 and 3.4. Fig. 4(a) displays the GB energy as a function of the misorientation angle between two grains (θ) for the lowest-energy atomic structures. The GBs examined do not exhibit a strong dependence of GB energy on θ. In particular, the GBs in the [001] system have similar GB energies in the range of 0.31–0.38 J m⁻² for 16.26° ≤ θ ≤ 77.32°. The Σ3(111) GB exhibits an exceptionally low GB energy of 0.07 J m⁻², as it is a twin boundary with a near-bulk atomic environment.

Fig. 4 GB properties as functions of the misorientation angle θ: (a) the GB energy, (b) the excess volume and (c) the coordination-number change per GB area. The quantities at 300 K were derived from the average atomic positions described in the caption of Fig. 3.

Fig. 4(b) shows the excess volume as a function of θ. For the [001] system, the excess volume varies slightly with θ, with a cusp for θ = 61.9°. The [11̄0] system exhibits a complex dependence of excess volume on θ. The GBs for θ > 70.53° tend to have larger excess volumes than those for θ < 70.53°, while their GB energies remain similar except for the Σ3(111) GB. The correlation coefficients between the GB energies and excess volume are 0.38 and 0.46 for the [001] and [11̄0] systems, respectively, indicating only a weak correlation; excess volume is therefore not a suitable metric for predicting the energetics of PbTe GBs. A similar conclusion was reported for fcc metals.⁷¹

The excess volumes at 300 K are smaller than those at 0 K for all GBs examined. If the GB regions expand with the same extent as the bulk, the excess volumes at the two temperatures should be identical. Thus, this result indicates that the thermal expansion of the GB regions is smaller than that in the bulk. It is also noteworthy that several GBs exhibit notable decreases in excess volume from 0 K to 300 K. This primarily results from changes in atomic structures, as shown in Fig. 3.

Fig. 4(c) shows the coordination-number change as a function of θ. It also varies with θ and is correlated with the excess volume, with correlation coefficients of 0.81 and 0.76 for the [001] and [11̄0] systems, respectively. The coordination-number change exhibits notable differences between 0 K and 300 K for several GBs, as is the case for the excess volume. These differences also arise from differences between their atomic structures at the two temperatures.

3.2 GB thermal conductivity and its anisotropy

This section investigates the structural dependence of κ_l for the GBs listed in Table 1. We evaluated κ_l along three axial directions as shown in Fig. 1(c). Fig. 5(a) shows the dependence of κ_l on θ for the [001] system. κ_l along the GB normal (z-axis), which is related to the GB thermal resistance, exhibited the most significant reduction as expected. The reduction from the bulk value depends on the GB structure and ranges from 68% to 81% (0.71 and 0.42 W m⁻¹ K⁻¹). κ_l for the [11̄0] system are also plotted in Fig. 5(b). Similarly, κ_l is lowest along the direction normal to the GB plane, with a reduction ranging from 39% to 80% (1.37 and 0.44 W m⁻¹ K⁻¹). The Σ3(111) GB (with θ = 70.53°) exhibits a notably higher κ_l owing to its high symmetry; however, even in this case, κ_l is reduced by approximately 40% in all three directions. When κ_l normal to the GB plane are converted to thermal resistances R_GB according to eqn (5), the values range from 1.4 to 9.7 m² K G⁻¹ W⁻¹, corresponding to a difference of approximately a factor of seven (Fig. 5(c)). This clearly demonstrates that lattice thermal conduction is strongly dependent on the GB structure.

Fig. 5 (a) and (b) Lattice thermal conductivity κ_l at 300 K for the GB models with rotation axes of (a) [001] and (b) [11̄0]. κ_l along three axial directions, namely normal to the GB planes (z), parallel to the rotation axis (x), and the other in-plane direction (y), are shown. The horizontal dashed line indicates κ_l for the perfect crystal. (c) GB thermal resistance estimated from κ_l along the GB-normal direction.

Qin et al. also evaluated κ_l of PbTe twin boundaries using an MLIP and the Green–Kubo approach.³⁷ The estimated R_GB is 1.1 m² K G⁻¹ W⁻¹, which is close to our value of 1.4. This suggests that, if the developed MLIPs are sufficiently accurate, the calculated R_GB would converge toward well-defined values. Issota et al. experimentally reported an average R_GB of 75 m² K G⁻¹ W⁻¹ in SnTe.²⁴ Our calculated R_GB values are about one order of magnitude smaller. Possible reasons include the fact that we consider relatively simple symmetric tilt GBs rather than complex random GBs and neglect effects such as GB segregation, in addition to the differences in κ and GB structures between SnTe and PbTe. A similar discrepancy has also been observed for Si GBs,^35,72 suggesting the necessity of treating more complex GBs in future studies.

In the directions parallel to the GB plane, κ_l decreased by 20–30% along the rotation axis (x-axis), and by 40–62% along the other in-plane direction (y-axis) for the [001] system (Fig. 5(a)). For the [11̄0] system, κ_l is reduced by 39–56% and 37–67% along the x- and y-axes, respectively (Fig. 5(b)). These large reductions show that phonon scattering along the GB plane is also non-negligible at nanometer GB spacing (e.g., nanocrystalline samples), although this effect is often neglected.

A major difference between the [001] and [11̄0] systems is that κ_l along the rotation axis is more significantly reduced for the [11̄0] system, indicating that thermal transport along [11̄0] is more sensitive to the GB effect. To elucidate the origin of this difference, the contributions of longitudinal and transverse waves on κ_l along the [001] and [11̄0] directions in the perfect crystal and GBs were also calculated, as shown in Fig. 6. In the perfect crystal (Fig. 6(a)), although total κ_l is the same for both directions, the longitudinal and transverse contributions differ significantly: along the [001] direction, longitudinal waves contribute to κ_l across the entire frequency range while transverse waves contribute only via low-frequency acoustic modes; along [11̄0] both longitudinal and transverse waves contribute across the spectrum. Accordingly, the longitudinal and transverse contributions account for 81% and 19% for [001], and 45% and 55% for [11̄0]. In the [001] system, the overall reduction in κ_l is relatively small because the dominant longitudinal component is less strongly scattered (Fig. 6b). In contrast, along [11̄0], both the longitudinal and transverse components decreased by comparable amounts. The [001] direction is parallel to the Pb–Te bonds, probably leading to the larger longitudinal contribution and the weaker influence of GB scattering. These results suggest that the influence of the GB on lattice thermal conduction along the planes depends significantly on the crystallographic orientation.

Fig. 6 (a) Spectral κ_l along [001] and [11̄0] directions in the perfect crystal at 300 K. The contributions from longitudinal and transverse waves on the spectral κ_l are also shown. (b) Contributions from longitudinal and transverse waves on κ_l for the GBs with [001] and [11̄0] rotation axes. The horizontal red and blue lines indicate the wave contributions for the perfect crystal.

3.3 Governing factors of GB thermal conductivity

We examined the factors governing κ_l, including GB energy, excess volume, coordination number, bond length, and other structural metrics. The averaged GB structures obtained from MD simulations at 300 K were used for the analysis (note that the correlation becomes less apparent when the 0 K structures are used for the analysis, Fig. S2).

Fig. 7(a) shows that the coordination-number change is well correlated with κ_l for along the GB-normal direction (see Fig. S3 for the result with an alternative definition of the coordination number). κ_l exhibits an exponential decay with increasing coordination-number change, with a large Pearson correlation coefficient r of −0.82. The [001] and [11̄0] GBs follow the same trend, indicating that the GB structure primarily determines κ_l in this direction. The sharp decrease in κ_l from the perfect crystal to the twin GB is likely attributable to grain misorientation rather than to the GB structure itself. A comparable correlation is found for the GB excess volume (Fig. S4), which demonstrates that reduced coordination and non-optimum packing of atoms at the GB core suppress phonon transport.

Fig. 7 Lattice thermal conductivity κ_l as a function of (a) coordination-number change and (b) bond-length distortion for the three axial directions.

The decrease in κ_l with increasing coordination-number change has also been reported for ionic MgO,^29,30 which adopts the same rocksalt structure as PbTe. On the other hand, for covalent Si, no clear correlation with excess volume has been found; instead, the GB energy has a moderate correlation with κ_l.²⁸ PbTe is considered to exhibit metavalent bonding, and reduced coordination at GBs may induce a more covalent character.²⁵ The similarity of the κ_l suppression mechanism to MgO implies that, for GB phonon scattering, the crystal structure and the resulting GB structures are more important than the bond character in the vicinity of the interface.

Along the x-axis (the rotation axis), excluding the twin boundary, κ_l tends to increase as the coordination-number change becomes larger (Fig. 7(a)). There is no physically reasonable mechanism by which the coordination-number change would lead to higher κ_l, and this apparent correlation should be indirect and mediated by other structural parameters. However, they have not been identified in the present analysis.

For the GB-parallel directions, bond-length distortions, defined as the sum of distortion indices divided by the GB area, were found to exhibit modest correlations with κ_l (Fig. 7(b)). Along the x-axis, the correlation depends on the rotation axis, with r = −0.87 for the [001] system and r = −0.54 for the [11̄0] system. Along y-axis, a common trend appears irrespective of the rotation axis with r = −0.70. These results indicate that, for the GB-parallel directions, κ_l is influenced by local structural modifications induced by GBs as well as crystallographic orientation. GB energy shows a similar distribution with the bond-length distortion, although the correlation with κ_l is weak probably because it is evaluated at 0 K (Fig. S2(d)).

To investigate the microscopic origins of the reduction and anisotropy in κ_l, we mapped the κ^atom_l,i on the zy plane. The results for the Σ5(210)/[001] GB in the three axial directions are shown in Fig. 8(a) and (b). These maps were obtained by smearing κ^atom_l,i with two different widths to facilitate interpretation of the atomically-resolved thermal conductivity and its spatial dependence. Along the GB-normal direction, κ^atom_l,i exhibited no dependence on the local structure, and all atoms showed uniformly low conductivity. This indicates that the GB acts as a rate-limiting barrier for the overall heat flux, resulting in a reduced effective thermal conductivity even within the grains. This behavior is consistent with the observation in Fig. 7(a) that the coordination-number change at the GBs governs κ_l. A similar tendency has also reported for ionic MgO^29,30 and SrTiO₃.³³ On the other hand, along the directions parallel to the GB plane, κ^atom_l,i clearly depends on the local structure: it is smaller near the GB and larger inside the grains. Because the GB structures are periodically repeated along the heat flux, κ^atom_l,i is strongly suppressed near the GB, whereas it becomes higher in the grain interiors where the crystal lattice is repeated.

Fig. 8 (a) and (b) Local thermal conductivities κ^local_l along the three axial directions in the Σ(210)/[001] GB, mapped with Gaussian broadening widths of (a) σ = 0.6 Å and (b) σ = 2.0 Å. (c) Atomic thermal conductivities κ^atom_l,i for all GB models as a function of the minimum distance from the two GB planes. Brighter colors indicate higher point densities.

Fig. 8(c) shows κ^atom_l,i in all GB models as a function of the distance from the GB plane for the three axial directions. The uniformly low κ^atom_l,i along the GB-normal direction is commonly observed, indicating that the GBs govern thermal transport in this direction over a range significantly longer than 5 nm. Along the GB-parallel directions, κ^atom_l,i increases with increasing distance from the GB, suggesting that local structural factors, such as elastic strains induced by the GB, suppress thermal transport over a few nm of the boundary. As κ^atom_l,i within the grains does not reach the bulk value for any of the GBs (including the twin), the GB influence probably persists up to a distance of approximately 10 nm along the GB-parallel directions. κ_l along the rotation axis (x-axis) is higher in the [001] system than in the [11̄0] system (Fig. 5(a) and (b)), and this trend is also observed in κ^atom_l,i profiles (Fig. 8(c)). These results again indicate that the spatial extent of GB influence on thermal transport depends on the crystallographic direction.

3.4 Relationship between local atomic environments and κ^atom_l,i

This section examines the distributions of local atomic environments (LAEs) embedded in the GBs and their correlations with κ^atom_l,i. The distributions of the coordination number and distortion index are shown in Fig. 9(a). For the [001] system, the coordination number lies between 4.5 and 6 and the distortion index is at most 0.08, whereas for the [11̄0] system these ranges broaden to 3.5–6.5 and up to 0.15, respectively. These comparisons suggest that [11̄0] GBs include more distorted LAEs than [001] GBs. Moreover, within the range of GBs analyzed in this study, the region covered by the [11̄0] GBs nearly coincides with that for all GBs examined, covering most LAEs seen in the [001] GBs. To confirm this, we also employed the SOAP descriptor, a more flexible structure descriptor that also incorporates angular information. To assess short- and long-range structural differences, we set the cutoff radii to 4.58 Å (at a distance of the second nearest neighbor) and 25 Å (approximately half the separation between the two GB planes), respectively, and the results confirm the same trend both cutoffs (Fig. S5). Such diverse distortions in the [11̄0] GBs may contribute to the lower κ_l along the GB-parallel directions, compared with the [001] system, as shown in Fig. 5(a) and (b).

Fig. 9 (a) Two-dimensional histogram of local atomic environments (LAEs) binned by coordination number and distortion index. Panels show the results for the [001] GBs, [11̄0] GBs, and all GBs combined. Color indicates LAE counts per bin on a logarithmic scale. (b) and (c) Atomic thermal conductivities κ^atom_l,i for all GB models plotted against (b) coordination number and distortion index and (c) local distortion factor with a SOAP cutoff of 25 Å. Brighter colors indicate higher point densities.

Before analyzing the structure–κ^atom_l,i relationship, the difference between the Pb and Te was evaluated. For the present NNP, κ^atom_l,i of Pb and Te in the perfect crystal are nearly identical to κ_l, 2.2 W m⁻¹ K⁻¹ at 300 K. In contrast, in the GB models the κ^atom_l,i of Te are larger than that of Pb by 0.13–0.20 W m⁻¹ K⁻¹ on average in each axial direction (Fig. S6). This difference can be attributed to the phonon frequencies of Pb and Te. As the atomic masses of Pb and Te are 207.2 and 127.6, respectively, lighter Te vibrates at higher frequencies than Pb. Correspondingly, spectral κ_l shows that Te makes a larger contribution in the higher-frequency range (Fig. S7). κ^atom_l,i for Te becomes larger probably because higher-frequency optical phonons are less scattered by GBs. The frequency-dependent scattering has also been reported in our previous studies for PbTe point defects³⁶ and MgO GBs.²⁹

Because the elemental dependence is minor compared with the effects of GB structure and crystallographic orientation, we hereafter discuss the structure–κ^atom_l,i relationship without separating Pb and Te. Fig. 9(b) shows that the correlation between the coordination number and κ^atom_l,i is negligible (r = 0.19), although the coordination-number change correlates with κ_l along the GB-normal direction (Fig. 7(a)). This is because the effective thermal conduction is suppressed even inside the grains by the long-range influence of GBs for this direction, as expected from Fig. 8. In contrast, the distortion index correlates modestly with κ^atom_l,i along the GB-parallel directions, with r ∼ −0.5 to −0.6 (Fig. 9(b)). κ^atom_l,i tends to decrease with increasing distortion index, as expected. However, atoms in the grain interior, where the distortion is small (<0.025), exhibit a wide range of κ^atom_l,i at comparable distortion levels. This indicates that short-range descriptors are insufficient to explain κ^atom_l,i even for the GB-parallel directions.

To evaluate short-range and long-range structural distortions on κ^atom_l,i, we defined a SOAP-based local distortion factor (LDF), as the Euclidean norm of the difference between SOAP vectors in the GB model and the perfect crystal. With a cutoff of 4.58 Å, the correlations with κ^atom_l,i are similar to those for distortion index, because this cutoff primarily probes the structural distortions of first neighbors (Fig. S8). By increasing the cutoff to 25 Å, long-range structural distortions induced by GBs are better captured, as the distributions that were concentrated near zero for distortion index and short-cutoff SOAP become more dispersed (Fig. 9(c)). However, the correlations are modestly improved only for y-axis (r = −0.67), whereas they remain low along the rotation axis (r = −0.55/−0.53 for [001]/[11̄0] GBs), and are negligible along the GB normal. The limited correlation even incorporating long-range structural distortions indicates the need for more macroscopic structural descriptors.

Previous MLIP studies on point defects and solid solutions in the PbTe system^36,73 suggest that variations in atomic force constants affect κ_l. Our recent work on point defects³⁶ quantified the difference in the harmonic phonon density of states (DOS) between atoms near a defect and in the bulk. This metric reflects changes in the local force field beyond the nearest neighbors and is directly related to the phonons that transport heat. However, its correlation with κ^atom_l,i was weak along the three axes for all GBs, similar to the distortion index (Fig. S9). This again highlights the necessity of macroscopic structural descriptors, particularly for extended defects with higher dimensionality than point defects.

From these observations, we finally consider how to construct descriptors for predicting κ^atom_l,i and hence κ_l for GBs and microstructures. First, local, short-range structure around a given atom is insufficient for prediction, even with large cutoffs of a few nm. This is because phonon mean free paths can be long (up to micrometers) and defects influence thermal transport over extended distances. Therefore, in addition to microscopic descriptors, macroscopic descriptors must incorporate the nano- and mesoscale distribution and network of defects. For GBs in particular, the spatial range over which a GB effectively suppresses thermal transport differs between the GB-normal and GB-parallel directions and further depends on the crystallographic orientation of grains. Thus, beyond quantifying macroscopic disorder, it is essential to evaluate an anisotropic, heat-flux-direction-dependent effective structural disorder. If such microscopic and macroscopic descriptors are evaluated, and large dataset of κ^atom_l,i are compiled for various GBs and polycrystalline models, machine-learning approaches such as graph neural networks should enable prediction of κ^atom_l,i and, by averaging them, the macroscopic κ_l. This framework would support design guidelines for tailoring thermal transport through defect engineering.

4 Conclusions

In this work, we systematically calculated κ_l and its anisotropy for 24 symmetric tilt PbTe GBs with [001] and [11̄0] rotation axes based on NNP-based PMD simulations. Along the GB-normal direction, κ_l varies depending on the coordination-number change occurring near the GBs, because the GBs strongly limit thermal transport even within the grains. In contrast, along the GB-parallel directions, κ_l tends to vary not with the coordination-number change but with the bond-length distortion and the crystallographic orientation of grains. The spatial extent of the GB effect likely reaches several hundred nanometers for the GB-normal direction, depending on the phonon mean free path, whereas that for the GB-parallel direction is limited to about 10 nm. κ^atom_l,i shows almost no correlation with local structure along the GB-normal direction, and even along the GB-parallel directions its behavior cannot be sufficiently explained in terms of distortion index and other local descriptors. These findings suggest that more macroscopic descriptors, such as the structures, spatial distribution, and network of GBs with respect to the direction of heat flux, are required for predicting κ_l. Combining high-accuracy MLIP-based datasets of κ_l and κ^atom_l,i for diverse GBs and polycrystals with more sophisticated methods for analyzing both local and macroscopic structures will pave the way for exploiting the GB structures and morphologies to further improve the thermoelectric performance of PbTe.

Author contributions

Susumu Fujii: conceptualization, data curation, funding acquisition, investigation, resources, software, validation, visualization, writing – original draft, writing – review & editing. Tatsuya Yokoi: conceptualization, data curation, funding acquisition, investigation, resources, software, validation, visualization, writing – original draft, writing – review & editing. Katsuyuki Matsunaga: resources, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). The code for NNP-PMD simulations can be found at https://github.com/NU-programs/NNP. The version of the NNP employed for this study is version 7. Supplementary information: supplementary tables, fgures, and GB structures in the POSCAR format. See DOI: https://doi.org/10.1039/d5ta09725h.

Acknowledgements

This work was supported by JST FOREST (Grant Number JPMJFR235X and JPMJFR2464) and JSPS KAKENHI (Grant Number JP22H04914, JP23K13544, JP23K04381, JP23H01671, and JP23K26365). The computations in this study were partly carried out using the resource offered by Research Institute for Information Technology, Kyushu University (under the category of General Projects) and by the Supercomputer Center in the Institute for Solid State Physics, the University of Tokyo.

References

1. G. J. Snyder and E. S. Toberer, Nat. Mater., 2008, 7, 105–114.
2. J. He and T. M. Tritt, Science, 2017, 357, eaak9997.
3. K. Biswas, J. He, I. D. Blum, C. I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid and M. G. Kanatzidis, Nature, 2012, 489, 414–418.
4. Y. Xiao and L.-D. Zhao, npj Quantum Mater., 2018, 3, 55.
5. Y. Zhong, J. Tang, H. Liu, Z. Chen, L. Lin, D. Ren, B. Liu and R. Ang, ACS Appl. Mater. Interfaces, 2020, 12, 49323–49334.
6. B. Jia, D. Wu, L. Xie, W. Wang, T. Yu, S. Li, Y. Wang, Y. Xu, B. Jiang, Z. Chen, Y. Weng and J. He, Science, 2024, 384, 81–86.
7. M. H. Lee, J. H. Yun, G. Kim, J. E. Lee, S.-D. Park, H. Reith, G. Schierning, K. Nielsch, W. Ko, A.-P. Li and J.-S. Rhyee, ACS Nano, 2019, 13, 3806–3815.
8. A. A. El-Sharkawy, A. M. Abou El-Azm, M. I. Kenawy, A. S. Hillal and H. M. Abu-Basha, Int. J. Thermophys., 1983, 4, 261–269.
9. G. Akhmedova and D. S. Abdinov, Inorg. Mater., 2009, 45, 854–858.
10. Y. Pei, X. Shi, A. Lalonde, H. Wang, L. Chen and G. J. Snyder, Nature, 2011, 473, 66–69.
11. T. Shiga, J. Shiomi, J. Ma, O. Delaire, T. Radzynski, A. Lusakowski, K. Esfarjani and G. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 155203.
12. C. W. Li, O. Hellman, J. Ma, A. F. May, H. B. Cao, X. Chen, A. D. Christianson, G. Ehlers, D. J. Singh, B. C. Sales and O. Delaire, Phys. Rev. Lett., 2014, 112, 175501.
13. Y. Lu, T. Sun and D.-B. Zhang, Phys. Rev. B, 2018, 97, 174304.
14. J. He, J. R. Sootsman, S. N. Girard, J.-C. Zheng, J. Wen, Y. Zhu, M. G. Kanatzidis and V. P. Dravid, J. Am. Chem. Soc., 2010, 132, 8669–8675.
15. B. Xiang, J. Liu, J. Yan, M. Xia, Q. Zhang, L. Chen, J. Li, X. Y. Tan, Q. Yan and Y. Wu, J. Mater. Chem. A, 2019, 7, 18458–18467.
16. B. Jia, Y. Huang, Y. Wang, Y. Zhou, X. Zhao, S. Ning, X. Xu, P. Lin, Z. Chen, B. Jiang and J. He, Energy Environ. Sci., 2022, 15, 1920–1929.
17. S.-H. Lo, J. He, K. Biswas, M. G. Kanatzidis and V. P. Dravid, Adv. Funct. Mater., 2012, 22, 5175–5184.
18. P. Jood, M. Ohta, A. Yamamoto and M. G. Kanatzidis, Joule, 2018, 2, 1339–1355.
19. Y. Li, D. Mei, H. Wang, Z. Yao, T. Zhu and S. Chen, Mater. Lett., 2015, 140, 103–106.
20. P. Xu, W. Zhao, X. Liu, B. Jia, J. He, L. Fu and B. Xu, Adv. Mater., 2022, 34, 2202949.
21. Z. Chen, Z. Jian, W. Li, Y. Chang, B. Ge, R. Hanus, J. Yang, Y. Chen, M. Huang, G. J. Snyder and Y. Pei, Adv. Mater., 2017, 29, 1606768.
22. Y. Wu, P. Nan, Z. Chen, Z. Zeng, R. Liu, H. Dong, L. Xie, Y. Xiao, Z. Chen, H. Gu, W. Li, Y. Chen, B. Ge and Y. Pei, Adv. Sci., 2020, 7, 1902628.
23. J. J. Kuo, S. D. Kang, K. Imasato, H. Tamaki, S. Ohno, T. Kanno and G. J. Snyder, Energy Environ. Sci., 2018, 11, 429–434.
24. E. Isotta, S. Jiang, G. Moller, A. Zevalkink, G. J. Snyder and O. Balogun, Adv. Mater., 2023, 35, 2302777.
25. R. Wu, Y. Yu, S. Jia, C. Zhou, O. Cojocaru-Mirédin and M. Wuttig, Nat. Commun., 2023, 14, 719.
26. H. Zhang, M. Shen, C. Stenz, C. Teichrib, R. Wu, L. Schäfer, N. Lin, Y. Zhou, C. Zhou, O. Cojocaru-Mirédin, M. Wuttig and Y. Yu, Small Sci., 2025, 5, 2300299.
27. C.-L. Fu, M. Cheng, N. T. Hung, E. Rha, Z. Chen, R. Okabe, D. C. Carrizales, M. Mandal, Y. Cheng and M. Li, Adv. Mater., 2025, 37, 2505642.
28. J. Hickman and Y. Mishin, Phys. Rev. Mater., 2020, 4, 033405.
29. S. Fujii, T. Yokoi and M. Yoshiya, Acta Mater., 2019, 171, 154–162.
30. S. Fujii, T. Yokoi, C. A. Fisher, H. Moriwake and M. Yoshiya, Nat. Commun., 2020, 11, 1854.
31. S. Fujii, K. Funai, T. Yokoi and M. Yoshiya, Appl. Phys. Lett., 2021, 119, 231604.
32. Y. Liu, Y. Bian, A. Chernatynskiy and Z. Han, Int. J. Heat Mass Tran., 2019, 145, 118791.
33. S. Fujii, H. Isobe, W. Sekimoto and M. Yoshiya, Scr. Mater., 2025, 258, 116524.
34. P. Friederich, F. Häse, J. Proppe and A. Aspuru-Guzik, Nat. Mater., 2021, 20, 750–761.
35. S. Fujii and A. Seko, Comput. Mater. Sci., 2022, 204, 111137.
36. T. Yokoi, S. Fujii, Y. Ogura and K. Matsunaga, Phys. Rev. B, 2025, 112, 024111.
37. M. Qin, X. Zhang, J. Zhu, Y. Yang, Z. Ti, Y. Shen, X. Wang, X. Liu and Y. Zhang, J. Mater. Chem. A, 2023, 11, 10612–10627.
38. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
39. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
40. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
41. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
42. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
43. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
44. M. Parrinello and A. Rahman, Phys. Rev. Lett., 1980, 45, 1196–1199.
45. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
46. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
47. J. Behler, Int. J. Quantum Chem., 2015, 115, 1032–1050.
48. J. Behler, Chem. Rev., 2021, 121, 10037–10072.
49. N. Artrith, A. Urban and G. Ceder, Phys. Rev. B, 2017, 96, 014112.
50. S. Shah, F. Palmieri and M. Datum, Neural Netw., 1992, 5, 779–787.
51. T. B. Blank and S. D. Brown, J. Chemom., 1994, 8, 391–407.
52. T. Yokoi, Y. Kondo, K. Ikawa, A. Nakamura and K. Matsunaga, J. Mater. Sci., 2021, 56, 3183–3196.
53. T. Yokoi, H. Kato, Y. Oshima and K. Matsunaga, J. Phys. Chem. Solids, 2023, 173, 111114.
54. M. Yoshiya, A. Harada, M. Takeuchi, K. Matsunaga and H. Matsubara, Mol. Simul., 2004, 30, 953–961.
55. J. H. Irving and J. G. Kirkwood, J. Chem. Phys., 1950, 18, 817.
56. Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio and A. Harju, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 094301.
57. A. J. Gabourie, Z. Fan, T. Ala-Nissila and E. Pop, Phys. Rev. B, 2021, 103, 205421.
58. S. Fujii, M. Yoshiya and C. A. J. Fisher, Sci. Rep., 2018, 8, 11152.
59. S. Fujii, M. Yoshiya, A. Yumura, Y. Miyauchi, M. Tada and H. Yasuda, J. Electron. Mater., 2014, 43, 1905–1915.
60. S. Fujii and M. Yoshiya, J. Electron. Mater., 2016, 45, 1217–1226.
61. H.-S. Yang, G.-R. Bai, L. Thompson and J. Eastman, Acta Mater., 2002, 50, 2309–2317.
62. L.-F. Huang, B. Grabowski, E. McEniry, D. R. Trinkle and J. Neugebauer, Phys. Status Solidi B, 2015, 252, 1907–1924.
63. M. Wagih and C. A. Schuh, Scr. Mater., 2023, 237, 115716.
64. W. H. Baur, Acta Crystallogr., Sect. B, 1974, 30, 1195–1215.
65. L. Himanen, M. O. J. Jäger, E. V. Morooka, F. Federici Canova, Y. S. Ranawat, D. Z. Gao, P. Rinke and A. S. Foster, Comput. Phys. Commun., 2020, 247, 106949.
66. J. Laakso, L. Himanen, H. Homm, E. V. Morooka, M. O. Jäger, M. Todorović and P. Rinke, J. Chem. Phys., 2023, 158, 234802.
67. M. A. Caro, Phys. Rev. B, 2019, 100, 024112.
68. J. P. Darby, J. R. Kermode and G. Csányi, npj Comput. Mater., 2022, 8, 166.
69. Y. Zhou, J.-Y. Yang, L. Cheng and M. Hu, Phys. Rev. B, 2018, 97, 085304.
70. M. Huang, P. Zhai, S. I. Morozov, W. A. Goddard III, G. Li and Q. Zhang, J. Alloys Compd., 2023, 959, 170429.
71. D. L. Olmsted, S. M. Foiles and E. A. Holm, Acta Mater., 2009, 57, 3694–3703.
72. T. Harada, K. Kutsukake, N. Usami, T. Ikari and A. Fukuyama, J. Appl. Phys., 2024, 136, 205703.
73. K. Conley, C. Gerber, A. Novick, T. Berriodi, E. S. Toberer and A. J. Karttunen, Mater. Horiz., 2025, 12, 8084–8094.

---