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

It is well-known that the calculation of thermal conductivity using classical molecular dynamics (MD) simulations strongly depends on the choice of the appropriate interatomic potentials. As proven for the case of graphene, while most of the available interatomic potentials estimate the structural and elastic constants with high accuracy, when employed to predict the lattice thermal conductivity they however lead to a variation of predictions by one order of magnitude. Here we present our results on using machine-learning interatomic potentials (MLIPs) passively fitted to computationally inexpensive ab-initio molecular dynamics trajectories without any tuning or optimizing of hyperparameters. These first-attempt potentials could reproduce the phononic properties of different two-dimensional (2D) materials obtained using density functional theory (DFT) simulations. To illustrate the efficiency of the trained MLIPs, we consider polyaniline C3N nanosheets. C3N monolayer was selected because the classical MD and different first-principles results contradict each other, resulting in a scientific dilemma. It is shown that the predicted thermal conductivity of 418 ± 20 W mK−1 for C3N monolayer by the non-equilibrium MD simulations on the basis of a first-attempt MLIP evidences an improved accuracy when compared with the commonly employed MD models. Moreover, MLIP-based prediction can be considered as a solution to the debated reports in the literature. This study highlights that passively fitted MLIPs can be effectively employed as versatile and efficient tools to obtain accurate estimations of thermal conductivities of complex materials using classical MD simulations. In response to remarkable growth of 2D materials family, the devised modeling methodology could play a fundamental role to predict the thermal conductivity.


Introduction
Thermal conductivity is among the most important intrinsic properties of materials playing a key role in energy device engineering and thermal management [1][2][3][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][6][7][8][9][10][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 −1 [15][16][17][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 −1 [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][29][30][31][32][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][35][36][37][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 3 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.

Computational methods
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 3 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.

Ab-initio modelling
First-principles calculations in this work were conducted using the Vienna Ab-initio Simulation Package (VASP) [44][45][46]. For the considered 2D systems, we used generalized gradient approximation (GGA) with Perdew -Burke − Ernzerhof (PBE) [47] method. For the bulk silicon and boron arsenide (BAs) we nonetheless found that ordinary GGA/PBE does not reproduce the experimental lattice constants very accurately and such that for those systems we employed the Perdew-Burke-Ernzerhof revised for solids (PBEsol) [7]. Plane-wave cutoff energies of 600 eV and 400 eV were considered for carbon-based and the rest of the systems, respectively. For geometry optimization, the convergence criterion for the energy and forces were set to (10 -5 ) eV and 0.001 eV/Å, respectively. Density functional perturbation theory (DFPT) simulations were performed over super-cell structures using VASP package. We then employed PHONOPY code [48] with the DFPT results as inputs to acquire phonon dispersions and group velocities as well. Ab-initio molecular dynamics (AIMD) simulations were performed with a time step of 1 fs, over super-cells suing the 3 × 3 × 1 and 2 × 2 × 2 Monkhorst-Pack [49] k-point gird for the monolayer and bulk systems, respectively. The AIMD simulation outputs include the atomic position and corresponding forces and total energy and stresses at every simulation time step, which were used to create the training data.

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: are the energy, atomic forces and stresses in the training set, respectively, and and σ MTP k,ij 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].

Molecular dynamics modelling
Non-equilibrium molecular dynamics (NEMD) simulations were carried out to predict the lattice thermal conductivity of C 3 N monolayer. We used Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [53] package along with the fitted MTP to introduce the atomic interactions. The details of NEMD calculations in this work are the same as those we have used in our previous studies [26,54,55]. In this approach periodic boundary conditions were applied along the planar directions, in which we used a small time increment of 0.25 fs. To simulate the steady-state heat transfer, we first relax the structures at the room temperature using the Nosé-Hoover thermostat method (NVT). Then few rows of atoms at the two ends were fixed and the rest of simulation box was divided into 22 slabs. Next a 20 K temperature difference was applied between the first (hot) and last (cold) slabs. In this process, the desired temperatures at the two ends were secured by the NVT method, while the remaining of the system was simulated using the constant energy method. The energy added (to hot slab) or removed (from cold slab) within the system by the NVT thermostat was used to calculate the steady state heat flux, which results in the formation of a temperature gradient. As in our previous works [26,54,55], lattice thermal conductivity was then calculated according to the applied heat flux and established temperature gradient using the one-dimensional form of the Fourier law.

Results and discussions
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 2 N [58] (figure 1(d)), C 3 N [59] or C 3 B (figure 1(f)) and C 7 N 6 [60] (figure 1(f)) and carbon-free MoS 2 ( figure 1(g)) and penta-SiP 2 ( figure 1(h)). In this study, we also considered the bulk boron arsenide ( figure 1(k)) and silicon ( figure 1(m)) lattices to .6392 Å and for SiP 2 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 2 monolayer. Phagraphene and Haeckelite exhibit more complicated lattices which can be representative of highly defective graphene sheets. C 3 N and C 3 B nanosheets are also symmetrical but they include two atomic types and thus they show more complicated dynamical behavior than graphene. C 7 N 6 nanosheet is a nanoporous, low symmetry binary system and therefore is the most complicated system considered in this work.
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 7 N 6 in which we considered a 2 × 3 × 1 lattice. For graphene, phagraphene, haeckelite, C 2 N, C 3 N, C 3 B, MoS 2 and SiP 2 , '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 2 N and SiP 2 ), 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 3 N and C 3 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.
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 2 , C 7 N 6 , C 2 N, SiP 2 (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 3 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 3 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][68][69][70].
The NEMD predictions for the length effect on the thermal conductivity of C 3 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 3 N, the thermal conductivity at infinite length would be measured to be 418 ± 20 W mK −1 . 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 3 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 −1 . In a latest NEMD study [79] authors also reported consistent values for the length effect on the thermal conductivity of C 3 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 −1 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 −1 for C 3 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 3 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 3 N thermal conductivity was reported via iteratively solving the BTE using the ShengBTE [29] package. Notably, the thermal conductivity of C 3 N by first principles were reported to be 380 [80], 128 [81], 80 [82], 380 [83] and 482 [84] W mK −1 , 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.
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 2 [85] nanomembranes, boron-nitride multilayers [86] or carbon-based 2D nanostructures [87][88][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 3 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 7 N 6 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.

Conclusion
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.