Efficient machine-learning based interatomic potentialsfor exploring thermal conductivity in two-dimensional materials

Thermal conductivity is among the most important intrinsic properties of materials playing a key role in energy device engineering and thermal management [1–4]. For the majority of advanced applications, as those in electronics and energy storage/conversion systems, components with higher thermal conductivities are highly desirable to enhance the heat dissipation rates and suppress the overheating issues. In this regard, materials with higher thermal conductivities can minimize the need for complicated and costly thermal management systems and may substantially improve efficiency. On the other side, for thermally insulating and thermoelectric materials, a component with the minimal thermal conductivity is more suitable to reduce the energy losses or to improve the thermal to electric energy conversion efficiency, respectively. Accordingly the accurate evaluation of a material's thermal conductivity is a fundamental issue for benchmarking relevant materials for a given application purpose [5–11]. While for the bulk materials there exist various experimental ways for estimating the thermal conductivity, when the size of the materials decreases and approach the nanoscale, this can turn into a challenging and complicated task.

Since the first experimental isolation of graphene [12, 13] in 2004, two-dimensional (2D) materials have been among the most attractive and vibrant class of materials, with wide-range of application prospects, from nanoelectronics to structural components. In particular thermal management and heat spreaders or harvesters are the current subject of great concern and efforts [1]. For the 2D materials with thicknesses in atomic scale, the experimental measurement of the thermal conductivity has been found to be quite a challenging task. For example, Sadeghi et al [14] reported several drawbacks of different experimental methods for measuring the thermal conductivity of 2D materials. In order to complement and validate the experimental characterization of the 2D materials thermal conductivity, accurate theoretical approaches are therefore paramount. In this regard, classical molecular dynamics (MD) simulations and density functional theory plus Boltzmann transport equation (DFT-BTE) are the most commonly employed methods to calculate the thermal conductivity of 2D materials. Nonetheless, the accuracy of MD estimations strongly depends on the accuracy of the interatomic potentials, and estimated thermal conductivities may substantially vary upon the choice of potentials. For instance, in the case of a monolayer graphene, experimental measurements obtain values within 1500–5300 W mK⁻¹ [15–18], while using the MD simulations on the basis of original Terosff [19], AIREBO [20], REBO [21] and optimized Tersoff [22], theoretical thermal conductivity (assuming graphene with 'infinite length') is found to be 870 [23], 709.2 [24], 350 [25] and ∼3000 W mK⁻¹ [26, 27], respectively, showing an order of magnitude variation. Nonetheless when using the first-principles DFT-BTE methods, the estimation for the thermal conductivity may also change depending on the choices for exchange-correlation functional and computational details for the 2nd and 3rd order force-constants calculations [8, 28–33].

In recent years, the astonishing advances in machine-learning techniques have opened new possibilities to address critical challenges in various fields, also in materials science [9, 34–38]. For example, actively trained machine-learning interatomic potentials [39] have been successfully employed to predict novel materials [40, 41] and examine lattice dynamics [42] and thermal conductivity [43] of bulk materials. As shown in the latest studies [42, 43], such machine learning-based interatomic potentials bring a DFT level accuracy for the computed energies and forces using classical molecular dynamics simulations. Nonetheless, the drawback of the current methodology is that it requires a rather complicated setup, in which a few hundred single-point DFT calculations have to be performed during the active learning training procedure.

Here we propose a simple, fast and accurate methodology on the basis of passively trained machine-learning interatomic potentials to estimate the thermal conductivity of 2D materials and related heterostructures. This process only requires inexpensive ab-initio molecular dynamics trajectories of less than 4 ps, without any need in additional DFT calculations. We trained interatomic potentials for several 2D materials, and those potentials are shown to reproduce very closely the phonon dispersions and phonon group velocities obtained by computationally costly DFT calculations. We then employed a fitted MLIP to estimate the thermal conductivity of polyaniline C₃N monolayer using the popular non-equilibrium molecular dynamics simulations. We argue that the proposed approach not only improves the accuracy of classical molecular dynamics simulations but may yield more reproducible estimates as compared with those by the first-principles DFT-BTE method. The prospects for the improvements of the proposed approach will be also briefly discussed.

As outlined in the introduction section, in this work we develop MLIPs passively trained over short AIMD trajectories. A developed MLIP was then used to practically predict the thermal conductivity of C₃N monolayer, by the classical MD simulations. To explain the employed methods, we respectively provide the details of first-principles calculations, development of MLIPs and non-equilibrium molecular dynamics simulations.

2.1. Ab-initio modelling

2.2. Training of interatomic potentials

We employed moment tensor potential (MTP) [40, 50] as an accurate and computationally efficient model of describing interatomic interaction. MTPs are based on the representation of atomic environments in the form of inertia tensors of various ranks multiplied by radial polynomial functions (Chebyshev polynomials). Such a representation allows for the implementation of many-body interactions. The functional form of the MTPs depends on the positions and types of neighboring atoms within the cutoff radius. This functional form is invariant with respect to orthogonal Euclidean transformations (rotations, reflections), as well as permutations of atoms of the same type. MTPs are parameterized by a set of coefficients that are found in the training procedure from the requirement to minimize the difference between the predicted and first principles energies, forces and stresses on a certain set of configurations (also called a training set), via:

$$\begin{equation}\mathop \sum \limits_{k = 1}^K \left[ {{w_e}{{\left( {E_k^{{\text{AIMD}}} - E_k^{{\text{MTP}}}} \right)}^2} + {w_f}\mathop \sum \limits_i^N \left| {f_{k,i}^{{\text{AIMD}}} - f_{k,i}^{{\text{MTP}}}} \right|{{\text{ }}^2} + {w_s}\mathop \sum \limits_{i,j = 1}^3 \left| {\sigma _{k,ij}^{{\text{AIMD}}} - \sigma _{k,ij}^{{\text{MTP}}}} \right|{{\text{ }}^2}} \right] \to {\text{min}},{\text{ }}\end{equation} \tag{ 1 }$$

where $E_k^{{\text{AIMD}}}$, $f_{k,i}^{{\text{AIMD}}}{\text{ }}$and $\sigma _{k,ij}^{{\text{AIMD}}}{\text{ }}$are the energy, atomic forces and stresses in the training set, respectively, and $E_k^{{\text{MTP}}}$, $f_{k,i}^{{\text{MTP}}}{\text{ }}$and $\sigma _{k,ij}^{{\text{MTP}}}$ are the corresponding values calculated with the MTP, K is the number of the configurations in the training set, N is the number of atoms and w_e, w_f and w_s are the non-negative weights that express the importance of energies, forces and stresses in the optimization problem, respectively, which in our study were set to 1, 0.1 and 0.001, respectively. The training procedure is thus a nonlinear optimization which is solved using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. MTPs have been tested on a number of atomic modeling applications [51], where they demonstrate excellent characteristics for reproducing the properties of both single-component [41, 52] and multicomponent materials [40, 43].

2.3. Molecular dynamics modelling

The focus of the current work is on the development of MLIP for the modeling of the thermal conductivity of 2D materials using the commonly used classical MD modeling. Nonetheless, this approach is also valid to simulate the thermal conductivity of bulk lattices. In figure 1, the studied structures in this study are illustrated, which include full carbon 2D structures of graphene (figure 1(a)), phagraphene [56] (figure 1(b)) and haeckelite [57] (figure 1(c)) and carbon-based binary lattices of C₂N [58] (figure 1(d)), C₃N [59] or C₃B (figure 1(f)) and C₇N₆ [60] (figure 1(f)) and carbon-free MoS₂ (figure 1(g)) and penta-SiP₂ (figure 1(h)). In this study, we also considered the bulk boron arsenide (figure 1(k)) and silicon (figure 1(m)) lattices to examine the accuracy of MTP for reproducing the phonic properties of bulk crystals. The hexagonal lattice constants of graphene, haeckelite, C₂N, C₃N, C₃B, C₇N₆ and MoS₂ were found to be 2.4674, 7.0869, 8.3277, 4.8605, 5.1735, 6.7949 and 3.1837 Å, respectively. The lattice parameters of phagraphene are 8.1084 and 6.6392 Å and for SiP₂ monolayer the lattice length is 4.8376 Å [61]. As mentioned earlier, these lattice constants are on the basis of PBE and they may change slightly with different functional or plane-wave cutoff energy. For the bulk boron arsenide and silicon, we measured lattice constants of 4.779 and 5.436 Å using the PBEsol functional, respectively, which match closely with experimentally measured values of 4.777 [62, 63] and 5.43 Å, respectively. In the supporting information document, we include the unit-cells of energy minimized lattices for the convenience of future studies. Among the studied systems, graphene includes the most symmetrical lattice, followed by bulk silicon and boron arsenide and MoS₂ monolayer. Phagraphene and Haeckelite exhibit more complicated lattices which can be representative of highly defective graphene sheets. C₃N and C₃B nanosheets are also symmetrical but they include two atomic types and thus they show more complicated dynamical behavior than graphene. C₇N₆ nanosheet is a nanoporous, low symmetry binary system and therefore is the most complicated system considered in this work.

Zoom In Zoom Out Reset image size

The objective of this study is to develop passively learned MLIPs, which can be trained rapidly and later efficiently employed in the classical MD simulations. To this end, the training set should include possible representative configurations which may occur during the classical MD simulations in order to avoid extrapolation. On this basis, for the creation of training data sets, AIMD simulations were carried out at different temperatures from 50 to 800 K, with overall simulation times of less than 4 ps. The low-temperature simulations can be useful to capture low displacement modes, equivalent with long wave-length acoustic phonons, whereas the simulations at higher temperatures are useful to better describe high-frequency optical modes. In addition, since in the NEMD simulations relatively long free-standing 2D sheets will be included, the AIMD trajectories at higher temperatures are necessary to include large local deformations that may happen during the lengthy calculations. We also note that AIMD simulations were conducted over n × n × 1 super-cells for monolayer, with an exception for C₇N₆ in which we considered a 2 × 3 × 1 lattice. For graphene, phagraphene, haeckelite, C₂N, C₃N, C₃B, MoS₂ and SiP₂, 'n' was equal to 6, 2, 2, 3, 3, 3, 5 and 3, respectively, and for bulk BAs and Si rectangular super-cells with 144 atoms were considered. The inclusion of larger super-cells can be useful especially for 2D materials in order to better describe the local deflections. With larger super-cells the total time of required AIMD calculations for the training can be decreased. Using the AIMD results, we trained MTPs with 117 and 189 parameters for the monoelemental and binary systems, respectively. In figure 2 (find figure S1 (stacks.iop.org/JPMATER/3/02LT02/mmedia) for C₂N and SiP₂), we compared the phonon dispersion relations predicted by passively trained MTPs with those by the DFT. We note that cumulative correlation factor can be also used to evaluate the agreement between the first-attempt MTP and the DFT results [64]. As it is clear from the comparison, the agreement with DFT results is spectacular, particularly for the low-frequency acoustic modes. This is a highly significant outcome since our first-attempt MTPs were trained automatically and no further modifications/optimization was applied/performed. One notes however that for C₃N and C₃B, the optical modes with the highest frequencies are less accurately captured by MLIP. For those cases increasing the MTPs constants can actually improve the accuracy. For silicon for instance, we compare the results with those by minimal-Tersoff [65], which has been very recently parameterized and provides the most accurate results for the silicon thermal conductivity. It is conspicuous that the fitted MLIP outperforms the minimal Tersoff [65], when reproducing the phonon dispersion. It is a noticeable finding since basically developed MLIPs were trained only over AIMD trajectories.

Zoom In Zoom Out Reset image size

Besides the phonon dispersions, reproducing the phonons group velocity can be also an important indication in accurately estimating the lattice thermal conductivity by the classical MD simulations. In figure 3 we thus compare the phonon group velocities predicted by MTPs with those obtained with DFT. Remarkably, the passively trained MTPs can very closely reproduce the phonon group velocities for the studied samples, particularly for graphene, MoS₂, C₇N₆, C₂N, SiP₂ (see figure S2) and BAs. Notably, the phonon group velocities for the low-frequency acoustic modes have been very accurately reproduced, which is a highly promising observation, since these modes are usually the main heat carriers in 2D materials. It is expected that by including more parameters in the fitted MLIPs the agreement for the phonon group velocity will be greatly improved (find figure S3 for the case of silicon). For this aim other methodologies like the learning on-the-fly [39, 40] can be also considered. Comparing the phonon dispersion and group velocity results, it is apparent that C₃N monolayer exhibits the maximal disagreement in the MLIP and DFT estimations. This structure was selected for the evaluation of thermal conductivity by the NEMD simulations. We remind that the C₃N was fist synthesized by Mahmood et al [59] in 2016 and can be also considered among the most attractive carbon-based 2D materials, which has been extensively explored recently for various applications [67–70].

Zoom In Zoom Out Reset image size

The NEMD predictions for the length effect on the thermal conductivity of C₃N monolayer at room temperature on the basis of a passively learned MTP are plotted in figure 4. As the characteristic of 2D materials with high thermal conductivities, the thermal conductivity at small lengths initially increases sharply due to the ballistic thermal transport. Normally the increase in the thermal conductivity at higher lengths is reduced and finally converges to reach the diffusive heat transfer regime. As a common approach the thermal conductivity at infinite length, k_∞, can be calculated by the extrapolation of the NEMD results for the samples with finite lengths, k_L, using the first-order rational curve fitting via 1/k_L = (1 + Λ/L)/k_∞ [71, 72], where Λ is the effective phonon mean free path. By assuming a thickness of 3.2 Å [54] for the single-layer C₃N, the thermal conductivity at infinite length would be measured to be 418 ± 20 W mK⁻¹. Optimized Tersoff [22] potential parameterized by Lindsay and Broido is currently accepted as the most accurate interatomic potential to study the thermal conductivity of graphene and in-general 2D carbon structures. On this basis, to simulate the thermal conductivity of carbon-nitride 2D systems, commonly the researchers use the Tersoff parameters constructed over Optimized Tersoff [22] proposed by the Kinarci et al [73]. Using the aforementioned interatomic potential setup along with the NEMD approach, the thermal conductivity of C₃N monolayer with infinite length was predicted to be 805–520 [54], 810–826 [74], 775 [75], 806 [76], 800 [77] and 780 [78] W mK⁻¹. In a latest NEMD study [79] authors also reported consistent values for the length effect on the thermal conductivity of C₃N monolayer with previous studies [54, 74, 76]. Although aforementioned NEMD estimations are based on the extrapolation of the small length data points to the infinity, different reports are all around 800 W mK⁻¹ and very close, which pronounce the reproducibility of this approach. As an exception, Dong et al [76] used the ReaxFF potential and estimated a thermal conductivity of 462 W mK⁻¹ for C₃N monolayer, however the same interatomic potential could not reproduce accurately the graphene thermal conductivity. On the other side and to the best of our knowledge, the thermal conductivity of C₃N has been also calculated using the DFT-BTE approach by four different groups. In the those studies, GGA/PBE functional and VASP package were employed for the DFT calculations, and the final value of the monolayer C₃N thermal conductivity was reported via iteratively solving the BTE using the ShengBTE [29] package. Notably, the thermal conductivity of C₃N by first principles were reported to be 380 [80], 128 [81], 80 [82], 380 [83] and 482 [84] W mK⁻¹, which shows six-fold scattering. This is an unexpected observation, not only because the first-principles estimations are expected to be with high level accuracy but also due to the fact that the simulation packages of aforementioned studies were basically the same. Clearly on one side, the majority of NEMD studies using a widely trusted interatomic potential lead to a convincingly reproducible estimation, but they also considerably overestimate the thermal conductivity by first principles. On the other side, the striking disagreement between different first-principles reports makes the understanding of the nanomaterial thermal conductivity unknown. It can be thus concluded that NEMD estimation on the basis of a MTP could be considered as a promising estimation, not only because it benefits from the classical MD reproducibility but also due to DFT level accuracy in the energy and force calculations.

Zoom In Zoom Out Reset image size

We note that in comparison with classical potentials like Tersoff, the computational costs of MTPs are by between 1 to 2 orders of magnitude higher. The high computational costs of MTPs could be mitigated by using a graphical processing units implementation. Further studies are required to examine the efficient number of constants and optimal cutoff distance for the MTPs in the evaluation of thermal conductivity. These modifications may easily accelerate the computational speed by several folds. In principle, NEMD simulations are computationally expensive only for materials with high thermal conductivities, because of the fact that simulations have to be conducted over long samples prior to extrapolating to the infinite length. It is thus believed that MLIP with numbers of constants between 100 and 200 parameters are practically usable with limited computational resources to evaluate the thermal conductivity of the majority of systems. On the other side we should note that the developed potentials due to the passive fitting methodology over pristine and defect-free samples may not be transferable for similar compositions or other problems, such as simulating bond breakages during the large deformations. Nevertheless, MLIPs offer outstanding flexibilities and enable active learning and tuning or optimizing of hyperparameters to improve the accuracy and transferability, but our study highlights that even with fitting over short AIMD trajectories and including limited number of parameters, the developed potentials can yet outperform the widely trusted/employed classical interatomic potentials. Notably, MLIPs offer novel possibilities not only to precisely predict the phonon dispersions but also to reproduce more complicated phononic properties, like phonon group velocity. However, finding efficient passive fitting approaches require more elaborated studies. Such a methodology has a large versatility and could be employed to explore the optimization of thermal conductivities in two-dimensional materials based compounds such as MoS₂ [85] nanomembranes, boron-nitride multilayers [86] or carbon-based 2D nanostructures [87–89]. Worthy to remind that when using the DFT-BTE approach, by decreasing the symmetry in the atomic lattice and depending on the cutoff distance, the number of required single-point DFT calculations increase substantially [29]. For-example, to acquire the 3rd order force-constants for the graphene and C₃N with a cutoff distance to 11th neighbors, one needs to conduct ∼200 and 700 DFT force-calculations, respectively. For low-symmetrical, nanoporous and anisotropic lattices, like C₇N₆ and graphyne/graphdiyne [90] structures, which normally show low thermal conductivities, we do believe that the proposed methodology can exhibit a unique performance to predict the thermal conductivity conveniently and accurately, with the need for minimal computation resources.

In summary, we have shown that machine-learning interatomic potentials (MLIP) trained over short ab-initio molecular dynamics trajectories (AIMD) are able to reproduce the phonon dispersions and phonon group velocity of complex 2D materials in remarkable agreements with first-principles results. Our classical estimation on the basis of a first-parameterized MLIP for the thermal conductivity of a 2D structure, for which different first-principles and classical reports considerably contradict each other, highlights the practical usefulness of the proposed approach, which simultaneously benefits from the reproducibility of the classical molecular dynamics modeling as well as density functional theory level accuracy in the force and energy evaluation. Nonetheless, more elaborated studies are required to find efficient number of parameters in MLIP, cutoff distance and corresponding prescription for the time and temperature of the AIMD simulations. Because of the fact that AIMD simulations are available by vast number of free and open-source first-principles packages, it is believed that the passive fitting can be efficiently employed to rapidly develop accurate and computationally efficient MLIPs for a desirable structure and later use it for the evaluation of thermal conductivity. This confirms a unique step within the classical molecular dynamics modeling, because it enables researchers to conveniently develop accurate interatomic potentials suitable for the modeling of the lattice thermal conductivity of a specific material using the classical MD, which would otherwise require a rather complicated fitting procedure. Taking into account the astonishing growth of 2D materials, which includes diverse compositions, with anisotropic, low-symmetrical and nanoporous members, we do believe that proposed methodology can be employed as a convenient and efficient approach to estimate the phononic thermal conductivity.

BM and XZ appreciate the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy within the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453). EVP, ISN, and AVS were supported by the Russian Science Foundation (Grant No. 18-13-00479). ICN2 is supported by the Severo Ochoa program from Spanish MINECO (Grant No. SEV-2017-0706) and funded by the CERCA Programme/Generalitat de Catalunya.

The following are the supplementary data to this article: Additionally, details of passive training of MTPs, VASP input files, energy minimized structures in VASP POSCAR format are accessible via: http://dx.doi.org/10.17632/7ppcf7cs27.1.