Equivariant graph neural network interatomic potential for Green-Kubo
  thermal conductivity in phase change materials
Thermal conductivity is a fundamental material property that plays an
essential role in technology, but its accurate evaluation presents a challenge
for theory. In this work, we demonstrate the application of $E(3)$-equivariant
neutral network interatomic potentials within Green-Kubo formalism to determine
the lattice thermal conductivity in amorphous and crystalline materials. We
apply this method to study the thermal conductivity of germanium telluride
(GeTe) as a prototypical phase change material. A single deep learning
interatomic potential is able to describe the phase transitions between the
amorphous, rhombohedral and cubic phases, with critical temperatures in good
agreement with experiments. Furthermore, this approach accurately captures the
pronounced anharmonicity that is present in GeTe, enabling precise calculations
of the thermal conductivity. In contrast, the Boltzmann transport equation
including only three-phonon processes tends to overestimate the thermal
conductivity by approximately a factor of 2 in the crystalline phases.
Sung-Ho Lee, Jing Li, Valerio Olevano, Benoit Sklénard
# Equivariant graph neural network interatomic potential for

###### Abstract

Thermal conductivity is a fundamental material property that plays an essential role in technology, but its accurate evaluation presents a challenge for theory. In this letter, we demonstrate the application of E(3)-equivariant neutral network interatomic potentials within Green-Kubo formalism to determine the lattice thermal conductivity in amorphous and crystalline materials. We apply this method to study the thermal conductivity of germanium telluride (GeTe) as a prototypical phase change material. A single deep learning interatomic potential is able to describe the phase transitions between the amorphous, rhombohedral and cubic phases, with critical temperatures in good agreement with experiments. Furthermore, this approach accurately captures the pronounced anharmonicity present in GeTe, enabling precise calculations of thermal conductivity. In contrast, the Boltzmann transport equation tends to overestimate it by approximately a factor of two in the crystalline phases.

Thermal conductivity is an intrinsic material property with deep implications in technology since it determines thermal management in the design of electronic devices [1; 2], and specifies the figure of merit in thermoelectric devices [3; 4]. Lattice vibrations, i.e. phonons, dominate heat transport in semiconductors and insulators. Much effort has been devoted to accurate calculations of lattice thermal conductivities from a microscopic perspective. The Boltzmann transport equation (BTE) [5; 6; 7], non-equilibrium Green function (NEGF) theory [8; 9; 10], and the Green-Kubo formula (GK) [11; 12; 13] are the three major approaches to lattice thermal conductivity calculations. BTE evaluates the response of phonon occupation to a temperature gradient, typically including three-phonon scattering processes, which limits its application to weakly anharmonic crystalline materials. NEGF treats phonons quantum mechanically and takes into account contact-channel interface scatterings and phonon anharmonicity by self-energies. However, it is computationally expensive [8]. GK provides the lattice thermal conductivity from the heat flux in an equilibrium molecular dynamics (MD) simulation, accounting for anharmonic effects to all orders [14]. Furthermore, recent developments extend GK to low temperatures [12; 13], which makes it a robust approach for a wide range of temperatures and materials. GK theory provides a unified approach to compute the lattice thermal conductivity in ordered and disordered solids. For harmonic amorphous systems, thermal transport can be described by the Allen and Feldman (AF) theory [15]. However, it has been shown that AF theory may be inadequate when anharmonic effects become important [16; 17].

The MD simulation in the GK approach requires a relatively long simulation time (up to a few nanoseconds) for adequate statistical sampling and an accurate description of interactions among atoms. Such long simulation times are affordable for MD with empirical force fields, but at the price of reduced accuracy and universality. _Ab initio_ MD has better accuracy but is too computationally expensive for large systems or long MD simulations. Extrapolation schemes have been proposed [18] to reduce the computational cost, but they are unsuitable for disordered solids.

In recent years, machine learning (ML) has emerged as a viable alternative for tasks that _ab initio_ methods have faced challenges with. In particular, machine learning interatomic potentials (MLIP) have been successful in predicting energies, forces and stress tensors orders of magnitude faster than first-principle methods, while retaining their accuracy. Thermal transport GK calculations have been reported with MLIP's relying on descriptor-based approaches, such as Behler-Parrinello neural networks (NN) or kernel-based methods [19; 20; 21]. Graph NN (GNN) interatomic potentials based on message passing architectures (MPNN) [22; 23; 24; 25] have been proposed as an alternative to hand-crafted descriptors, whereby structures are encoded as a graph with atoms represented as nodes that are connected by edges. In initial models, the information at nodes and edges of the GNN was made _invariant_ with respect to the Euclidean group \(E(3)\) (i.e. the group of translations, rotations and inversions in Euclidean space), and the atomic representations were limited to scalar interatomic distances [22]. Such models have since been generally superseded by MPNN architectures built on convolution operations that are _equivariant_ with respect to the E(3) group. In equivariant approaches, isometric transformations on the relative atomic displacement vector inputs are propagated through the network to correspondingly transform the outputs. Equivariant approaches have been shown to achieve substantially improved data efficiency and un
precedented accuracy compared to their invariant counterparts [23; 24; 25]. In MPNNs, many-body interactions are captured by iteratively propagating information along the graph at each layer in the network. This has the effect of extending the local receptive field of an atom to significantly beyond the cutoff radius, which renders parallelization impractical [26]. Recently, a strictly local equivariant neural network approach has been proposed to address this drawback [26]. In this architecture, information is stored as a per-pair quantity, and instead of nodes exchanging information with its neighbours via edges, a convolution operation acts on the cutoff sphere in the form of a set of invariant (scalar) latent features and a set of equivariant (tensor) latent features that interact at each layer.

In this letter, we demonstrate that the strictly local E(3)-equivariant NN can be employed to compute the temperature-dependent thermal conductivity of germanium telluride (GeTe) in various phases using GK theory. GeTe is a chalcogenide material employed in many technological applications, such as phase change nonvolatile memory storage [27; 28], thermoelectricity [3; 29; 30] and spintronics [31; 32; 33]. It undergoes a ferroelectric phase transition from the low temperature rhombohedral \(\alpha\)-GeTe (spacegroup \(R3m\)) to a cubic \(\beta\)-GeTe (spacegroup \(Fm\bar{3}m\)) at a Curie temperature of \(T_{c}\approx 650-700\) K [34; 35; 36]. Amorphous GeTe also plays an important role in technological applications. Therefore, GeTe is an ideal prototype phase change material for the study of lattice thermal conductivity using GK theory.

The thermal conductivity tensor within GK theory is defined as:

\[\kappa_{\alpha\beta}(T)=\frac{1}{k_{\text{B}}T^{2}V}\lim_{\tau\to\infty}\int_ {0}^{\tau}dt\,\langle j_{\alpha}(t)\cdot j_{\beta}(0)\rangle_{T}, \tag{1}\]

where \(k_{\text{B}}\) is the Boltzmann constant, \(T\) the temperature, \(V\) the volume, \(j_{\alpha}(t)\) the \(\alpha\)-th Cartesian component of the macroscopic heat flux, and \(\langle j_{\alpha}(t)\cdot j_{\beta}(0)\rangle_{T}\) the heat flux autocorrelation function (HFACF), with the symbol \(\langle\cdot\rangle_{T}\) denoting ensemble average over time and over independent MD trajectories.

The total heat flux of a system of \(N\) atoms is defined

\[\mathbf{j}(t)=\sum_{i=1}^{N}\frac{d}{dt}\left(\mathbf{r_{i}}E_{i}\right), \tag{2}\]

where \(E_{i}=m_{i}\mathbf{v}_{i}^{2}/2+U_{i}\) is the total energy (i.e. kinetic and potential energy) of atom \(i\) with mass \(m_{i}\), velocity \(\mathbf{v}_{i}\) and atomic positions \(\mathbf{r_{i}}\). In MLIPs, the partitioning \(E=\sum_{i}E_{i}\) of the total energy of the system into atomic contributions \(E_{i}\) allows the total heat flux of a periodic system to be expressed as [11] :

\[\mathbf{j}(t)=\sum_{i=1}^{N}\mathbf{v}_{i}E_{i}-\sum_{i=1}^{N}\sum_{j\neq i}\mathbf{r}_{ ij}\left(\frac{\partial U_{i}}{\partial\mathbf{r}_{ij}}\cdot\mathbf{v}_{j}\right) \tag{3}\]

where the sum over \(j\) runs over the atoms that are within the cutoff radius \(r_{c}\) of atom \(i\) defined for the MLIP. We implemented the calculation of Eq. (3) in the LAMMPS code [37]. The term \(\partial U_{i}/\partial\mathbf{r}_{ij}\) is obtained by automatic differentiation of atomic energies \(U_{i}\) computed by the MLIP. It was also used for the calculation of the virial tensor [11; 38], which is required to perform simulations in the isothermal-isobaric (NpT) ensemble.

To generate the reference dataset to train the MLIP, _ab initio_ MD simulations based on density functional theory (DFT) were performed with temperatures ranging from 100 K to 2500 K using the VASP code [39; 40]. The generalized gradient approximation of Perdew-Burke-Ernzerhof (PBE) [41] was used for the exchange-correlation energy and Grimme's D3 dispersion correction [42] was applied. The supercells contained 192 and 216 atoms for the initial rhombohedral and cubic structures, respectively. Then, 6000 structures in total were taken from the MD trajectories and recomputed to obtain more accurate energy, forces, and stress tensors. We used an energy cutoff of 400 eV and a \(2\times 2\times 2\)\(k\)-mesh to sample the Brillouin zone. The equivariant NN model was trained on energy, forces and stress using the Allegro package [26]. The root mean squared errors (RMSE) and mean absolute errors (MAE) on the predicted energies, forces and stress tensors on the test dataset are 0.90 meV/atom, 29.87 meV/A, 0.28 meV/A\({}^{3}\) and 1.07 meV/atom, 42.97 meV/A, 0.37 meV/A\({}^{3}\), respectively (see Suppl. Mat. for more information on the training procedure and dataset partitioning).

To further validate the MLIP, the equilibrium geometries of crystalline GeTe were optimized using the MLIP. For \(\alpha\)-GeTe, the lattice parameter was \(a=4.42\) A and the angle \(\alpha=57.13^{\circ}\), close to DFT results of \(a=4.41\) A and \(\alpha=57.42^{\circ}\). Similarly, for \(\beta\)-GeTe, the MLIP yields \(a=4.24\) A, in excellent agreement with the lattice parameter from DFT of \(a=4.23\) A.

Moreover, the phonon dispersion from the MLIP is in excellent agreement with DFT for both \(\alpha\) and \(\beta\)-GeTe, as shown in Fig. 1. In particular, our model describes optical phonons well, which is usually challenging for MLIPs [19; 43]. Imaginary soft phonon modes in cubic GeTe are also well described by the MLIP, which is essential to capture the phase transition [44; 45]. These phonon dispersions were computed using the finite displacement method implemented in Phonopy [46] with \(3\times 3\times 3\) and \(5\times 5\times 2\) supercells of the conventional unit cells for cubic and rhombohedral phases, respectively. For the DFT calculations, we used the same settings as those used to generate the reference dataset. LO-TO splitting was not included in our calculations as long-range Coulomb interactions tend to be screened by free carriers in real samples [47].

We investigated the lattice dynamics of GeTe through MD simulations across the \(\alpha\to\beta\) phase transition with
our MLIP. For each temperature, GeTe supercells were first equilibrated for at least 200 ps in the NpT ensemble at ambient pressure with a 2 fs timestep in order to obtain the averaged temperature-dependent structural parameters shown in Fig. 2. The rhombohedral lattice parameter \(a\) and angle \(\alpha\) reach cubic values at \(T\approx 650\) K, in good agreement with experimental data.

By employing the temperature-dependent effective-potential (TDEP) method [48; 49; 50], the temperature-dependent interatomic force constants (IFCs) were extracted from a 600 ps MD simulation in the microcanonical ensemble, after equilibrating the system in the NVT ensemble using the structural parameters depicted in Fig. 2. By utilizing these IFCs, we computed phonon spectra as a function of temperature (refer to the Suppl. Mat. for more detailed information). Fig. 3 presents the evolution of the longitudinal and transverse optical phonon modes (\(\Gamma_{6}\) and \(\Gamma_{4}\), respectively) as a function of temperature. The softening of these two modes up to the Curie temperature is corroborated by previous theoretical studies [44; 45] and is comparable to experiments [51; 47; 52]. Beyond 650 K, the optical phonons merge, indicating the transition to the cubic phase where optical phonons exhibit three-fold degeneracy.

To compute the GK thermal conductivity of cubic, rhombohedral and amorphous GeTe, MD simulations with the MLIP were performed at different temperatures. The amorphous GeTe structure was generated using a melt-quench process (see Suppl. Mat.). The heat flux was calculated during MD simulations in the microcanonical ensemble and the ensemble average was performed over independent trajectories of at least 1 ns after equilibration in the NpT ensemble. After testing the convergence with respect to system size (see Suppl. Mat.), we used supercells containing 360 atoms for the rhombohedral phase and 512 atoms for the amorphous and cubic phases.

Figure 1: Comparison of phonon dispersions computed with DFT and with the MLIP of (a) \(\alpha\)-GeTe and (b) \(\beta\)-GeTe

Figure 3: Temperature evolution of A\({}_{1}\) and E optical phonon modes computed with the TDEP method and compared against experimental data from Ref. [51; 47; 52].

Figure 2: Evolution of (a) the lattice parameter \(a\) and (b) the angle \(\alpha\) as a function of temperature in the NpT MD simulations of crystalline GeTe, compared against experimental data from Ref. [34; 35; 36]. Simulated lattice parameters in (a) were shifted by \(-0.1\) Å.
Although cubic GeTe is metastable below \(T_{c}\), GK is able to determine its lattice thermal conductivity as it becomes dynamically stable at \(T\geq 300\) K (see finite temperature phonon spectra in Suppl. Mat.). Rhombohedral GeTe shows a higher thermal conductivity than cubic GeTe before 650 K (see Fig. 4) after which the two curves merge, reflecting the \(\alpha\rightarrow\beta\) phase transition.

The comparison against experiments is challenging because experimental values of lattice thermal conductivities of crystalline GeTe show a large dispersion. There are two reasons for this. First, thermal conductivity comprises a lattice contribution and an electronic contribution. Therefore, experimental lattice thermal conductivity is an indirect measurement, which is obtained by removing the electronic contribution, typically evaluated using the Wiedmann-Franz law that introduces an additional approximation from the Lorenz number. Second, the sample quality varies. Extrinsic scatterings due to defects may alter the thermal conductivity measurements. For example, an extra phonon-vacancy scattering has to be included in order to recover a good agreement with experimental data [54, 55]. Despite the significant experimental variations mentioned above, the calculated GK thermal conductivity values are found to fall within the range of experimental values.

The GK lattice thermal conductivity for the amorphous phase (solid green line) is in excellent agreement with the experimental data of Ref. [54] (green squares). This can be regarded as a direct comparison with the experiment since the electronic contribution to the thermal conductivity was found to be negligible in amorphous GeTe [56]. A previous study obtained a similar value of \(0.27\pm 0.05\) W\(\cdot\)m\({}^{-1}\cdot\)K\({}^{-1}\) at 300 K from GK simulations with a Behler-Parrinello-type MLIP [21]. The predicted thermal conductivity for amorphous GeTe is constant until \(\sim 450\) K. It then starts to increase, indicating a transition to a crystalline phase, as evidenced by the evolution of the radial distribution function (see Suppl. Mat.) and consistent with the amorphous-crystalline phase transition temperature observed experimentally [54].

To obtain the BTE thermal conductivity, we used the TDEP 2nd and 3rd order IFCs from MD simulations and a \(30\times 30\times 30\)\(q\)-mesh. This allows a direct comparison between GK and BTE as both calculations were on the same footing, with identical interatomic potential and the same temperature; the only difference being the thermal transport formalism. BTE overestimates the thermal conductivity by about 1.8 W\(\cdot\)m\({}^{-1}\cdot\)K\({}^{-1}\), which is about twice the GK result at 300 K, and about three times that at 900 K. Such overestimation is an indication that BTE cannot capture the strong anharmonicity exhibited by GeTe.

In conclusion, we developed an equivariant graph neural network interatomic potential to study thermal transport in amorphous and crystalline GeTe. The potential describes GeTe at a near-_ab initio_ level of accuracy for the rhombohedral, cubic and amorphous phases with a single model. Our potential also correctly captures phase transitions with Curie temperatures in good agreement with experimental data. Combined with the Green-Kubo theory, it can determine the lattice thermal conductivity not only for strongly anharmonic crystals, but also for the amorphous phase.

We thank F. Bottin and J. Bouchet for discussions about TDEP calculation. This work was performed using HPC/AI resources from GENCI-IDRIS (Grant 2022-A0110911995) and was partially funded by European commission through ECSEL-IA 101007321 project StorAIge and the French IPCEI program.
