LETTER • OPEN ACCESS You may also like - First-principles investigation of the valley Efficient machine-learning based interatomic and electrical properties of carbon-doped - graphyne-like BN sheet potentialsfor exploring thermal conductivity in two- Bo Chen, , Xiang-Qian Li et al. - Two-dimensional carbon nitride C6N dimensional materials nanosheet with egg-comb-like structure and electronic properties of a semimetal A Bafekry, M Shahrokhi, A Shafique et al. To cite this article: Bohayra Mortazavi et al 2020 J. Phys. Mater. 3 02LT02 - Hydrogen storage on Li-decorated B4N: a first-principle calculation insight Yuhua Wei, Feng Gao, Jiguang Du et al. View the article online for updates and enhancements. Recent citations - Thermal conductivity of straight-chain polytetrafluoroethylene: A molecular dynamics study Xinjian Liu et al - The MLIP package: moment tensor potentials with MPI and active learning Ivan S Novikov et al - Machine-learning interatomic potentials enable first-principles multiscale modeling of lattice thermal conductivity in graphene/borophene heterostructures Bohayra Mortazavi et al This content was downloaded from IP address 79.152.138.100 on 11/01/2022 at 11:10 J. Phys. Mater. 3 (2020) 02LT02 https://doi.org/10.1088/2515-7639/ab7cbb Journal of Physics: Materials LETTER Efficient machine-learning based interatomic potentialsfor OPEN ACCESS exploring thermal conductivity in two-dimensional materials RECEIVED 21 January 2020 Bohayra Mortazavi1, Evgeny V Podryabinkin2, Ivan S Novikov2, Stephan Roche3,4, Timon Rabczuk5, REVISED Xiaoying Zhuang1 and Alexander V Shapeev2 26 February 2020 1 ACCEPTED FOR PUBLICATION Institute of Continuum Mechanics, Leibniz Universität Hannover, Appelstraße 11, 30157, Hannover, Germany 2 4 March 2020 Skolkovo Institute of Science and Technology, Skolkovo Innovation Center, Nobel St. 3, Moscow 143026, Russia 3 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193, Barcelona, Spain PUBLISHED 4 9 April 2020 ICREA Institució Catalana de Recerca i Estudis Avancats, 08010, Barcelona, Spain 5 College of Civil Engineering, Department of Geotechnical Engineering, Tongji University, Shanghai, People’s Republic of China Original content from E-mail: stephan.roche@icn2.cat and zhuang@ikm.uni-hannover.de this work may be used under the terms of the Keywords: machine learning, thermal conductivity, molecular dynamics, density functional theory simulations, Two-dimensional Creative Commons polyaniline C3N monolayer Attribution 4.0 licence. Any further distribution Supplementary material for this article is available online of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Abstract 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 C3 N nanosheets. C3 N 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 C3 N 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. 1. Introduction 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 © 2020 The Author(s). Published by IOP Publishing Ltd J. Phys. Mater. 3 (2020) 02LT02 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−1 [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−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–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 C3 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. 2. 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 C3 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 First-principles calculations in this work were conducted using the Vienna Ab-initio Simulation Package (VASP) [44–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 2 J. Phys. Mater. 3 (2020) 02LT02 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. 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:   X K  X N X 3 we EAIMD k − Ek MTP 2 + wf fkAIMD ,i − fkMTP ,i 2 + ws σkAIMD ,ij − σkMTP ,ij 2 → min, (1) k=1 i i,j=1 where EAIMD k , fkAIMD ,i and σkAIMD ,ij are the energy, atomic forces and stresses in the training set, respectively, and MTP MTP MTP Ek , fk,i and σ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 we , wf and ws 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 Non-equilibrium molecular dynamics (NEMD) simulations were carried out to predict the lattice thermal conductivity of C3 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. 3. 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 C2 N [58] (figure 1(d)), C3 N [59] or C3 B (figure 1(f)) and C7 N6 [60] (figure 1(f)) and carbon-free MoS2 (figure 1(g)) and penta-SiP2 (figure 1(h)). In this study, we also considered the bulk boron arsenide (figure 1(k)) and silicon (figure 1(m)) lattices to 3 J. Phys. Mater. 3 (2020) 02LT02 Figure 1. Atomic structures of studied systems in this work. Top and side view of monolayer (a) graphene (b) phagraphene [56] (c) haeckelite [57] (d) C2 N [58] (e) C3 N [59] or C3 B (F) C7 N6 [60] (g) MoS2 (h) penta-SiP2 . 3D views of bulk (k) boron arsenide and (m) silicon. examine the accuracy of MTP for reproducing the phonic properties of bulk crystals. The hexagonal lattice constants of graphene, haeckelite, C2 N, C3 N, C3 B, C7 N6 and MoS2 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 SiP2 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 MoS2 monolayer. Phagraphene and Haeckelite exhibit more complicated lattices which can be representative of highly defective graphene sheets. C3 N and C3 B nanosheets are also symmetrical but they include two atomic types and thus they show more complicated dynamical behavior than graphene. C7 N6 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 C7 N6 in which we considered a 2 × 3 × 1 lattice. For graphene, phagraphene, haeckelite, C2 N, C3 N, C3 B, MoS2 and SiP2 , ‘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 C2 N and SiP2 ), 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 C3 N and C3 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 4 J. Phys. Mater. 3 (2020) 02LT02 Figure 2. Phonon dispersion relations of (a) graphene (b) phagraphene (c) haeckelite (d) C3 N (e) C3 B (F) C7 N6 (g) MoS2 and bulk (h) boron arsenide and (k) silicon. The continuous and black dotted lines represent the first-attempt MTP and DFT estimations, respectively. For the silicon (panel k) we also included the minimal-Tersoff [65] (black line) and experimental [66] (empty circles) results. 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, MoS2 , C7 N6 , C2 N, SiP2 (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 C3 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 C3 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]. The NEMD predictions for the length effect on the thermal conductivity of C3 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 5 J. Phys. Mater. 3 (2020) 02LT02 Figure 3. Phonon group velocity (v) of single-layer (a) graphene (b) phagraphene (c) haeckelite (d) C3 N (e) C3 B (F) C7 N6 (g) MoS2 (h) and bulk (h) boron arsenide and (k) silicon predicted by MTP and DFT. for the samples with finite lengths, kL , using the first-order rational curve fitting via 1/kL = (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 C3 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 C3 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 C3 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 C3 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 C3 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 C3 N thermal conductivity was reported via iteratively solving the BTE using the ShengBTE [29] package. Notably, the thermal conductivity of C3 N by first principles were reported to be 380 [80], 128 [81], 80 [82], 380 [83] and 482 [84] W mK−1 , which shows 6 J. Phys. Mater. 3 (2020) 02LT02 Figure 4. NEMD results (empty circles) for the length effect on the thermal conductivity of single-layer C3 N at 300 K. The line illustrates the best fit to the NEMD data points. 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 MoS2 [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 C3 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 C7 N6 and graphyne/graphdiyne [90] structures, which normally show low thermal 7 J. Phys. Mater. 3 (2020) 02LT02 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. 4. 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. Acknowledgments 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. Supplementary data 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. ORCID iDs Bohayra Mortazavi  https://orcid.org/0000-0003-3031-5057 Stephan Roche  https://orcid.org/0000-0003-0323-4665 Xiaoying Zhuang  https://orcid.org/0000-0001-6562-2618 References [1] Fu Y et al 2019 Graphene related materials for thermal management 2D Mater. 7 12001 [2] Renteria J, Nika D and Balandin A 2014 Graphene thermal properties: applications in thermal management and energy storage Appl. Sci. 4 525–47 [3] Sun S et al 2018 Improving thermal transport at carbon hybrid interfaces by covalent bonds Adv. Mater. Interfaces 5 1800318 [4] Nylander A N, Fu Y, Huang M and Liu J 2019 Covalent anchoring of carbon nanotube-based thermal interface materials using epoxy-silane monolayers IEEE Trans. Componen. Packag. Manuf. Technol. 9 427–33 [5] Balandin A A 2011 Thermal properties of graphene and nanostructured carbon materials Nat. Mater. 10 569 [6] Malekpour H et al 2014 Thermal conductivity of graphene laminate Nano Lett. 14 5155–61 [7] Csonka G I et al 2009 Assessing the performance of recent density functionals for bulk solids Phys. Rev. B 79 155107 [8] Huang S et al 2019 Anisotropic thermal conductivity in 2D tellurium 2D Mater. 7 15008 [9] Sun B and Barnard A S 2019 Visualising multi-dimensional structure/property relationships with machine learning J. Phys. Mater. 2 34003 [10] Chen Y, Ma J and Li W 2019 Understanding the thermal conductivity and Lorenz number in tungsten from first principles Phys. Rev. B 99 20305 8 J. Phys. Mater. 3 (2020) 02LT02 [11] Meng F, Ma J, He J and Li W 2019 Phonon-limited carrier mobility and temperature-dependent scattering mechanism of (C3 N)-SiC from first principles Phys. Rev. B 99 45201 [12] Novoselov K S et al 2004 Electric field effect in atomically thin carbon films Science 306 666–9 [13] Geim A K and Novoselov K S 2007 The rise of graphene Nat. Mater. 6 183–91 [14] Sadeghi M M, Pettes M T and Shi L 2012 Thermal transport in graphene Solid State Commun 152 1321–30 [15] Ghosh S et al 2008 Extremely high thermal conductivity of graphene: prospects for thermal management applications in nanoelectronic circuits Appl. Phys. Lett. 92 1–4 [16] Balandin A A et al 2008 Superior thermal conductivity of single-layer graphene Nano Lett. 8 902–7 [17] Jauregui L A et al 2010 Thermal transport in graphene nanostructures: experiments and simulations ECS 28 73–83 [18] Cai W et al 2010 Thermal transport in suspended and supported monolayer graphene grown by chemical vapor deposition Nano Lett. 10 1645–51 [19] Tersoff J 1989 Modeling solid-state chemistry: interatomic potentials for multicomponent systems Phys. Rev. B 39 5566–8 [20] Stuart S J, Tutein A B and Harrison J A 2000 A reactive potential for hydrocarbons with intermolecular interactions J. Chem. Phys. 112 6472–86 [21] Brenner D W et al 2002 A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons J. Phys. Condens. Matter. (https://doi.org/10.1088/0953-8984/14/4/312) [22] Lindsay L and Broido D A 2010 Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene Phys. Rev. B 81 205441 [23] Wei Z, Ni Z, Bi K, Chen M and Chen Y 2011 In-plane lattice thermal conductivities of multilayer graphene films Carbon NY 49 2653–8 [24] Hong Y, Ju M G, Zhang J and Zeng X C 2018 Phonon thermal transport in a graphene/MoSe2 van der Waals heterobilayer Phys. Chem. Chem. Phys. 20 2637–45 [25] Thomas J A, Iutzi R M and McGaughey A J H 2010 Thermal conductivity and phonon transport in empty and water-filled carbon nanotubes Phys. Rev. B 81 045413 [26] Mortazavi B and Rabczuk T 2015 Multiscale modeling of heat conduction in graphene laminates Carbon NY 85 1–7 [27] Fan Z et al 2017 Thermal conductivity decomposition in two-dimensional materials: application to graphene Phys. Rev. B 95 144309 [28] Lindsay L, Broido D A and Mingo N 2010 Flexural phonons and thermal transport in graphene Phys. Rev. B 82 115427 [29] Li W, Carrete J, Katcho N A, Mingo N and Sheng B T E 2014 A solver of the Boltzmann transport equation for phonons Comput. Phys. Commun. 185 1747–58 [30] Torres P, Alvarez F X, Cartoixà X and Rurali R 2019 Thermal conductivity and phonon hydrodynamics in transition metal dichalcogenides from first-principles 2D Mater. 6 35002 [31] Zeng Q et al 2019 Highly anisotropic thermoelectric properties of black phosphorus crystals 2D Mater. 6 45009 [32] Rashid Z, Nissimagoudar A S and Li W 2019 Phonon transport and thermoelectric properties of semiconducting Bi2Te2X (X = S, Se, Te) monolayers Phys. Chem. Chem. Phys. 21 5679–88 [33] Raeisi M, Ahmadi S and Rajabpour A 2019 Modulated thermal conductivity of 2D hexagonal boron arsenide: a strain engineering study Nanoscale 11 21799–810 [34] Oda H, Kiyohara S and Mizoguchi T 2019 Machine learning for structure determination and investigating the structure-property relationships of interfaces J. Phys. Mater. 2 34005 [35] Schleder G R, Padilha A C M, Acosta C M, Costa M and Fazzio A 2019 From DFT to machine learning: recent approaches to materials science–a review J. Phys. Mater. 2 32001 [36] Schmidt J, Marques M R G, Botti S and Marques M A L 2019 Recent advances and applications of machine learning in solid-state materials science Npj Comput. Mater. 5 83 [37] Dong Y et al 2019 Bandgap prediction by deep learning in configurationally hybridized graphene and boron nitride Npj Comput. Mater. 5 26 [38] Carleo G et al 2019 Machine learning and the physical sciences Rev. Mod. Phys. 91 45002 [39] Podryabinkin E V and Shapeev A V 2017 Active learning of linearly parametrized interatomic potentials Comput. Mater. Sci 140 171–80 [40] Gubaev K, Podryabinkin E V, Hart G L W and Shapeev A V 2019 Accelerating high-throughput searches for new alloys with active learning of interatomic potentials Comput. Mater. Sci 156 148–56 [41] Podryabinkin E V, Tikhonov E V, Shapeev A V and Oganov A R 2019 Accelerating crystal structure prediction by machine-learning interatomic potentials with active learning Phys. Rev. B 99 064114 [42] Ladygin V V, Korotaev P Y, Yanilkin A V and Shapeev A V 2020 Lattice dynamics simulation using machine learning interatomic potentials Comput. Mater. Sci 172 109333 [43] Korotaev P, Novoselov I, Yanilkin A and Shapeev A 2019 Accessing thermal conductivity of complex compounds by machine learning interatomic potentials Phys. Rev. B 100 144308 [44] Kresse G and Furthmüller J 1996 Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set Comput. Mater. Sci 6 15–50 [45] Kresse G and Furthmüller J 1996 Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set Phys. Rev. B 54 11169–86 [46] Kresse G 1999 From ultrasoft pseudopotentials to the projector augmented-wave method Phys. Rev. B 59 1758–75 [47] Perdew J, Burke K and Ernzerhof M 1996 Generalized gradient approximation made simple Phys. Rev. Lett. 77 3865–8 [48] Togo A and Tanaka I 2015 First principles phonon calculations in materials science Scr. Mater. 108 1–5 [49] Monkhorst H and Pack J 1976 Special points for Brillouin zone integrations Phys. Rev. B 13 5188–92 [50] Shapeev A V 2016 Moment tensor potentials: A class of systematically improvable interatomic potentials Multiscale Model. Simul. 14 1153–73 [51] Nyshadham C et al 2019 Machine-learned multi-system surrogate models for materials prediction Npj Comput. Mater. 5 51 [52] Novoselov I I, Yanilkin A V, Shapeev A V and Podryabinkin E V 2019 Moment tensor potentials as a promising tool to study diffusion processes Comput. Mater. Sci. 164 46–56 [53] Plimpton S 1995 Fast parallel algorithms for short-range molecular dynamics J. Comput. Phys. 117 1–19 [54] Mortazavi B 2017 Ultra high stiffness and thermal conductivity of graphene like C3N Carbon NY 118 25–34 [55] Mortazavi B et al 2019 Nanoporous graphene: A 2D semiconductor with anisotropic mechanical, optical and thermal conduction properties Carbon NY 147 377–84 9 J. Phys. Mater. 3 (2020) 02LT02 [56] Wang Z et al 2015 Phagraphene: A low-energy graphene allotrope composed of 5-6-7 carbon rings with distorted dirac cones Nano Lett. 15 6182–6 [57] Terrones H et al 2000 New metallic allotropes of planar and tubular carbon Phys. Rev. Lett. 84 1716–19 [58] Mahmood J et al 2015 Nitrogenated holey two-dimensional structures Nat. Commun. 6 1–7 [59] Mahmood J et al 2016 Two-dimensional polyaniline (C3N) from carbonized organic single crystals in solid state Proc. Natl Acad. Sci. 113 7414–19 [60] Mortazavi B, Shahrokhi M, Shapeev A V, Rabczuk T and Zhuang X 2019 Prediction of C7N6 and C9N4: stable and strong porous carbon-nitride nanosheets with attractive electronic and optical properties J. Mater. Chem. C 7 10908–17 [61] Shojaei F, Mortazavi B, Zhuang X and Azizi M 2020 Silicon diphosphide (SiP2) and silicon diarsenide (SiAs2): novel stable 2D semiconductors with high carrier mobilities, promising for water splitting photocatalysts Mater. Today Energy 16 100377 [62] Chu T L and Hyslop A E 1972 Crystal growth and properties of boron monoarsenide J. Appl. Phys. (https://doi.org/10.1063/1.1661106) [63] Buckeridge J and Scanlon D O 2019 Electronic band structure and optical properties of boron arsenide Phys. Rev. Mater. 3 51601 [64] Dong Y et al 2019 Molecular bridge thermal diode enabled by vibrational mismatch Phys. Rev. Appl. 11 24043 [65] Fan Z et al 2020 A minimal Tersoff potential for diamond silicon with improved descriptions of elastic and phonon transport properties J. Phys.: Condens. Matter. 32 135901 [66] Holt M et al 1999 Determination of phonon dispersions from X-Ray transmission scattering: the example of silicon Phys. Rev. Lett. 83 3317–19 [67] Shi L-B, Zhang -Y-Y, Xiu X-M and Dong H-K 2018 Structural characteristics and strain behaviors of two-dimensional C3N : first principles calculations Carbon NY 134 103–11 [68] Tagani M B and Vishkayi S I 2018 Polyaniline (C3N) nanoribbons: magnetic metal, semiconductor, and half-metal J. Appl. Phys. 124 84304 [69] Bafekry A, Stampfl C, Shayesteh S F and Peeters F M 2019 Exploiting the novel electronic and magnetic structure of C3N via functionalization and conformation Adv. Electron. Mater. 5 1900459 [70] Bafekry A, Stampfl C, Ghergherehchi M and Farjami Shayesteh S 2020 A first-principles study of the effects of atom impurities, defects, strain, electric field and layer thickness on the electronic and magnetic properties of the C2N nanosheet Carbon NY 157 371–84 [71] Schelling P K, Phillpot S R and Keblinski P 2002 Comparison of atomic-level simulation methods for computing thermal conductivity Phys. Rev. B 65 1–12 [72] Zhang X et al 2014 Thermal conductivity of silicene calculated using an optimized Stillinger-Weber potential Phys. Rev. B 89 054310 [73] KInacI A, Haskins J B, Sevik C and Ça�in T 2012 Thermal conductivity of BN-C nanostructures Phys. Rev. B 86 115410 [74] Hong Y, Zhang J and Zeng X C 2018 Monolayer and bilayer polyaniline C3N: two-dimensional semiconductors with high thermal conductivity Nanoscale 10 4301–10 [75] Cheng D H, Wang X, Ding W, Chen Y, Zhang J, Xin G and Cheng L 2018 Phonon thermal conduction in a graphene-C3N heterobilayer using molecular dynamics simulations Nanotechnology (https://doi.org/10.1088/1361-6528/aaf481) [76] Dong Y, Meng M, Groves M M, Zhang C and Lin J 2018 Thermal conductivities of two-dimensional graphitic carbon nitrides by molecule dynamics simulation Int. J. Heat Mass Transf. 123 738–46 [77] Song J et al 2019 Thermal conductivity of two-dimensional BC3: a comparative study with two-dimensional C3N Phys. Chem. Chem. Phys. (https://doi.org/10.1039/C9CP01943J) [78] Hatam-Lee S M, Rajabpour A and Volz S 2020 Thermal conductivity of graphene polymorphs and compounds: from C3N to graphdiyne lattices Carbon NY 161 816–26 [79] An M et al 2020 Mass difference and polarization lead to low thermal conductivity of graphene-like carbon nitride (C3N Carbon NY 162 202–8 [80] Wang H, Li Q, Pan H, Gao Y and Sun M 2019 Comparative investigation of the mechanical, electrical and thermal transport properties in graphene-like C3B and C3N J. Appl. Phys. 126 234302 [81] Kumar S, Sharma S, Babar V and Schwingenschlögl U 2017 Ultralow lattice thermal conductivity in monolayer C3N as compared to graphene J. Mater. Chem. A 5 20407–11 [82] Wang H et al 2018 Lone-pair electrons do not necessarily lead to low lattice thermal conductivity: an exception of two-dimensional penta-CN2 J. Phys. Chem. Lett. 9 2474–83 [83] Gao Y et al 2018 First-principles study of intrinsic phononic thermal transport in monolayer C3N Physica E 99 194–201 [84] Peng B et al 2018 Tuning thermal transport in C3N monolayers by adding and removing carbon atoms Phys. Rev. Appl. 10 34046 [85] Sledzinska M et al 2016 Thermal conductivity of MoS2polycrystalline nanomembranes 2D Mater. 3 035016 [86] Cai Q et al 2019 High thermal conductivity of high-quality monolayer boron nitride and its thermal expansion Sci. Adv. 5 eaav0129 [87] Bazrafshan S and Rajabpour A 2017 Thermal transport engineering in amorphous graphene: non-equilibrium molecular dynamics study Int. J. Heat Mass Transf. 112 379–86 [88] Rajabpour A, Bazrafshan S and Volz S 2019 Carbon-nitride 2D nanostructures: thermal conductivity and interfacial thermal conductance with the silica substrate Phys. Chem. Chem. Phys. 21 2507–12 [89] Tromer R M, Felix I M, Freitas A, Azevedo S and Pereira L F C 2019 Diboron-porphyrin monolayer: A new 2D semiconductor Comput. Mater. Sci 172 109338 [90] Mortazavi B et al 2019 N-, B-, P-, Al-, As-, Ga-graphdiyne/graphyne lattices: first-principles investigation of mechanical, optical and electronic properties J. Mat. Chem. C7 3025–36 10