This is the submitted version of the article: Mortazavi B., Podryabinkin E.V., Roche S., Rabczuk T., Zhuang X., Shapeev A.V.. Machine-learning interatomic potentials enable first-principles multiscale modeling of lattice thermal conductivity in graphene/borophene heterostructures. Materials Horizons, (2020). 7. : 2359 - . 10.1039/d0mh00787k. Available at: https://dx.doi.org/10.1039/d0mh00787k Machine Learning Interatomic Potentials Enable First-Principles Multiscale Modeling of Lattice Thermal Conductivity in Graphene/Borophene Heterostructures Bohayra Mortazavia,b, Evgeny V Podryabinkinc, Stephan Roche*d,f, Timon Rabczukg, Xiaoying Zhuanga,g** and Alexander V Shapeevc a Institute of Continuum Mechanics, Leibniz Universität Hannover, Appelstraße 11, 30157 Hannover, Germany. b Cluster of Excellence PhoenixD (Photonics, Optics, and Engineering–Innovation Across Disciplines), Gottfried Wilhelm Leibniz Universität Hannover, Hannover, Germany. c Skolkovo Institute of Science and Technology, Skolkovo Innovation Center, Nobel St. 3, Moscow 143026, Russia. d Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain. f ICREA Institució Catalana de Recerca i Estudis Avancats, 08010 Barcelona, Spain g College of Civil Engineering, Department of Geotechnical Engineering, Tongji University, Shanghai, China. Abstract One of the ultimate goals of computational modeling is to work accurately with the minimal needs of empirical information. First-principles approaches like density functional theory (DFT) simulations address the aforementioned goal, but they are limited to systems consisting of a few hundreds, or thousands of atoms at most. On the other hand, classical molecular dynamics (CMD) simulations and finite element method (FEM) are extensively employed to study larger and more realistic systems, but they are dependent on the empirical information Here, we highlight that machine-learning interatomic potentials (MLIPs) trained over short ab-initio molecular dynamics trajectories enable first-principles multiscale modeling, in which DFT simulations can be hierarchically bridged to efficiently simulate the macroscopic structures. As a case study, we explore the lattice thermal conductivity of a coplanar graphene/borophene heterostructures, recently synthesized experimentally (Sci. Adv. 2019; 5: eaax6444), for which no viable classical modeling alternative is presently available. Our MLIPs approach can efficiently predict the thermal conductivity of graphene and borophene pristine phases, as well as thermal conductance of polycrystalline morphologies of mixed graphene/borophene lattices. This study shows that MLIPs can be effectively and conveniently employed to achieve first-principles multiscale modeling via hierarchical employment of DFT/CMD/FEM simulations, which can play a critical role to design novel nanostructures. Corresponding authors: *stephan.roche@icn2.cat; **zhuang@ikm.uni-hannover.de From the engineering point of view, numerical simulations are currently an inevitable part of the design, which not only substantially reduce the final costs of a product but also enable 1 the optimization toward the enhancement of its functional performance. In comparison with conventional materials, experimental characterization techniques for the evaluation of nanomaterials and nanostructure properties are substantially more complicated, time- consuming and expensive as well. Like other engineered products, for the practical application of nanomaterials in various technologies, developments of accurate modeling approaches are critical to facilitate the design and further optimizations. In recent years the state of the art theoretical simulations have played a major role in the astonishing advances in the field of materials science. In this regard, modeling enables researchers to examine the stability and explore the properties of novel materials and structures purely through computer simulations. Notably, first-principles simulations can be currently employed to find possible synthesis routes for the design of novel materials1–3. As a recent example, boron nanosheets with different atomic lattices were epitaxially grown over the silver surface 4,5, a fabrication process that was originally proposed by the density functional theory (DFT) simulations6,7. The main drawback of first-principles DFT calculations is nonetheless related to their demanding computational cost, which limits the maximal size of studied systems to only several hundreds, or at most a few thousand atoms. Moreover, the computational costs of common DFT simulations normally scale exponentially by increasing the number of atoms, which rule out their practical application for studying large and disordered materials models. Classical molecular dynamics (CMD) simulation is among the most popular numerical approaches, and has been extensively employed to explore the properties and responses of nanostructured materials. Unlike the DFT simulations, the computational cost of CMD calculations scales linearly with the number of atoms so that the scale of several million atoms models can be reached. Nonetheless, the accuracy of CMD results strongly depends on the precision of the interatomic potentials in describing the energies and forces. As a well-known example, despite the bonding mechanism in the planar full sp2 carbon system, most of the currently available interatomic potentials cannot accurately reproduce the thermal conductivity of graphene. Moreover, for the novel materials and structures, it is a challenging task to find an interatomic potential that can access the stable structures, irrespective of the accuracy in estimating their basic properties. It is clear that in comparison with the DFT counterpart, the computational advantage of MD simulations comes with the costs of declined flexibility and accuracy as well. On the other side, continuum mechanic based method like the finite element method (FEM) offer robust solutions to study real engineering problems, but in these methods, the properties of the materials should be fully known prior to modeling design. It is thus conspicuous that for studying the properties and responses of nanomaterials, development of multiscale approaches, solving each method drawback, are crucially needed. The latest advances in the field of machine-learning methods have offered novel creative solutions to address critical challenges for diverse problems, especially in materials science 8– 12. In particular machine-learning interatomic potentials (MLIPs) have been recently successfully employed to predict novel materials 13,14 and examine lattice dynamics 15,16 and 2 thermal conductivity17 of conventional bulk materials. As proven in numerous recent studies15,17,18, MLIPs enable CMD simulations to be conducted within the DFT level accuracy for the computed energies and forces, but with computational costs scaling linearly with the number of atoms. Another critical advantage of MLIPs is that since they are elaborated from DFT simulations, they can be trained for a specific material composition and thus do not face the flexibility concern of standard CMD method. Accordingly MLIPs can provide unprecedented possibility to marry first-principles accuracy with multiscale modeling. To illustrate such strategy, we here examine the lattice thermal conductivity of graphene/borophene heterostructures19, as a truly challenging system to simulate accurately with conventional methods. To date, no classical interatomic potential that can accurately reproduce the structural properties of borophene nanosheets are available in the litterature. Moreover, worthy to remind that for a well-studied system like graphene, while the majority of interatomic potentials estimate the structural and elastic constants with good accuracy, when employed to predict the lattice thermal conductivity they usually lead to one order of magnitude variation. For the case of graphene, the experimental measured thermal conductivities are within 1500-5300 W/mK20–23, while CMD based estimates by the original Terosff 24, AIREBO 25, REBO 26 and optimized Tersoff 27, were reported to be 87028, 70929, 35030 and ~3000 W/m.K31,32, respectively. An accurate prediction of lattice thermal conductivity, using classical interatomic potentials, is therefore clearly a highly challenging issue. In addition, when using the CMD simulations to evaluate the interfacial thermal conductance, the interatomic potential must exhibit outstanding stability and accuracy as well, otherwise the calculations fail to simulate the steady-state heat transfer33–35. This issue is certainly very critical for the modeling of graphene/borophene interfaces, not only because of their different lattice structures but also due to the possibility of the formation of diverse types of grain boundaries. The main steps within the first-principles hierarchical multiscale modeling framework proposed here are summarized in Fig. 1. This includes four key steps: (1) DFT simulations; (2) development of MLIPs; (3) CMD simulations and (4) FEM modeling of effective lattice thermal conductivity. Within the DFT step, we first conducted the energy minimization of graphene and borophene lattices. Next, ten different possible grain boundaries (GBs) that may form between the graphene and borophene lattices (find Fig. 2a) were constructed. Since we conducted the DFT simulations within the plane-wave approach, the constructed models are periodic in all directions. Such that it is possible to construct two different grain boundaries in every DFT interface model. In order to create the required training sets for the development of MLIPs, ab-initio molecular dynamics (AIMD) simulations were conducted. These simulations were carried out for pristine phases (pure graphene or borophene) and heterostructures with geometry optimized interfaces at different temperatures of 100, 300, 600 and 700 K, each one with less than 1000 time steps. Since the AIMD trajectories are correlated within short time periods, subsamples of original trajectories were included in the training sets. In the next step, moment tensor potentials (MTPs)36 were parameterized to describe the interatomic interactions. Likely to classical counterparts, MTPs also include 3 parameters that are optimized over the training configurations provided by the AIMD simulations. In this work, two types of MTPs were developed, monoelemental potentials to simulate the pristine graphene or borophene and binary potentials to study the heterostructures. In the latter case, the created training sets not only include the AIMD trajectories from the interfaces but also those from the pristine phases. We first trained MTPs over subsampled AIMD trajectories. After the initial training of MTPs, the accuracy of the trained potentials were evaluated over the full AIMD trajectories and the configurations with high extrapolations grades18 were selected. The selected configurations were accordingly added to the original training sets and the final MTPs were developed by retraining of new clean potentials over the updated training sets. After the MTPs were trained, in the third step, they were employed to evaluate the thermal conductivity of pristine phases or calculate the interfacial thermal conductance of the grain boundaries via CMD simulations. In the last step, the effective lattice thermal conductivity was evaluated with the FEM method, in which the data provided by the third step were used to define the materials and interaction properties as well. Fig. 1, Main steps of first-principles multiscale modeling framework to simulate the lattice thermal conductivity of graphene/borophene heterostructures. In Fig. 2a, the atomic configurations of ten different graphene/borophene grain boundaries constructed in this study are illustrated. We remind that from a theoretical point of view, the grain boundaries in graphene are made from a series of pentagon/heptagon pairs37. Nonetheless for the case of MoS2 which also shows the hexagonal unit cell but with more complex geometry, according to high-resolution electron microscopy results, grain boundaries contain diverse forms of; pentagon-heptagon (5-7), tetragon-tetragon (4-4), tetragon-hexagon (4-6), tetragon-octagon (4-8) and hexagon-octagon (6-8) rings38–40. 4 Because of the different atomic structures of borophene and graphene and depending on the various tilting angles, it is expected that the graphene/borophene grain boundaries also show diverse configurations. From the constructed grain boundaries shown in Fig. 2a it is clear that they mainly include tetragon, pentagon, hexagon and octagon dislocations, but nonagon rings may also form as it occurred in the case of GB-1. We then conducted non-equilibrium molecular dynamics (NEMD) simulations to predict the interfacial thermal conductance and lattice thermal conductivity of pristine borophene. To this end we used the LAMMPS 41 package along with the trained MTPs to introduce the atomic interactions. In the NEMD approach periodic boundary conditions were applied along the planar directions using a simulation time step of 0.5 fs. As shown schematically in Fig. 2b, to simulate the steady-state heat transfer, we first relax the structures at room temperature using the Nosé-Hoover thermostat (NVT) method. Then few rows of atoms at the two ends were fixed and the rest of simulation box was divided into 22 slabs. Next a temperature difference of 20 K was applied between the first (hot) and last (cold) slabs. In this process, the desired temperatures at the two ends were controlled by the NVT method, while the remaining of the system was simulated without applying a thermostat. As shown in Fig. 2c for a sample of grain boundaries, to keep the applied temperature difference at every simulation time step an amount of energy is added to the hot slab and another amount of energy is removed from the cold slab by the NVT thermostat. As can be seen from Fig. 2c, the amounts of the energy added and removed to the system are accurately the same (that show linear patterns), confirming that the system has been precisely stayed under steady-state heat transfer condition. The slope of these energy curves can be used to calculate the applied steady-state heat flux (Hf). As shown in Fig. 2d, due to the existence of grain boundary, the temperature profile exhibits a sudden change at the interface (ΔT). It is noticeable that temperature gradient within the graphene region is negligible as compared with the borophene section, suggesting a considerably higher lattice thermal conductivity of graphene. The grain thermal conductance can be calculated as Hf/ΔT. For the pristine borophene, the temperature profile however illustrates a constant gradient, which can be used to estimate the thermal conductivity using the one-dimensional form of the Fourier law. In Fig. 2e, the calculated interfacial thermal conductance for the considered grain boundaries are compared. Notably, the thermal conductances of different grain boundaries are close. We also examined the length dependence and found that it does not affect the estimated thermal conductance, in agreement with a recent study on polycrystalline MoS2 42 and graphene/h-BN heterostructures43. These observations reveal that the interfacial resistance mainly stems from the very different phononic characteristics of graphene and borophene, in contrast with those of polycrystalline sheets in which the misorientation angle of adjacent sheets and density and type of dislocations cores play the critical role 42,43. Since the thermal conductance does not show substantial dependency on the types of formed defects in graphene/borophene interfaces, it is thus expected that including more extensive grain boundary configurations will not lead to considerable changes in the estimated effective lattice thermal conductivity. 5 Fig. 2, (a) Atomic configurations of constructed graphene/borophene grain boundaries (GB), (b) schematic illustration of non-equilibrium molecular dynamics (NEMD) method, (c) energy values added to the hot slab and removed from the cold slab by the NVT thermostat during every simulation time step, (d) established temperature profile showing a sudden drop at the interface (e) estimated interfacial thermal conductance of considered grain boundaries in panel (a). The length effect on the NEMD predictions for the lattice thermal conductivity of borophene monolayer at room temperature is plotted in Fig. 3a. Unlike the graphene, borophene shows anisotropic transport properties and such that the calculations were conducted along the armchair and zigzag directions. Sharp initial increases in the lattice thermal conductivity by increasing the sample length are observable. This length effect on the thermal conductivities suppresses at higher lengths and finally converges to reach the diffusive heat transfer regime. As a common approach, the thermal conductivity of borophene at infinite length, k∞, can be calculated by the extrapolation of the NEMD results for the samples with finite lengths, kL, using the first-order rational curve fitting via 1/kL=(1+Λ/L)/ k∞44,45, where Λ is the effective phonon mean free path. By assuming a thickness of 2.9 Å, the diffusive lattice thermal conductivity of borophene at room temperature along with the armchair and zigzag directions were estimated to be 52 and 112 W/mK, respectively. Another alternative to calculate the lattice thermal conductivity is to solve the Boltzmann transport equation. To that end we used the ShengBTE46 package, which offers full iterative solution of the Boltzmann transport equation to estimate the lattice thermal conductivity of graphene. The computationally demanding section of aforementioned approach is to acquire the third-order (anharmonic) interatomic force constants, which usually requires few hundred/thousand single-point DFT calculations over supercell lattices. In this work, second and third-order force constants were calculated using the density functional perturbation theory simulations and passively trained MTPs, respectively, over 10×10×1 super-cells (consisting of 200 atoms). For the evaluation of third-order anharmonic force constant, we considered the interactions with the eleventh nearest neighbours. In this case, by using the ShengBTE 46 package, we calculated the force constants using the MTP for 312 structures in a negligible time, which otherwise with DFT would require decent computational resources. On 6 the basis of passively trained MTP trained over AIMD simulations within the PBE/GGA functional, the diffusive lattice thermal conductivity of graphene was estimated to be 3600 W/mK, which falls within the experimentally measured values of 1500-5300 W/mK20–23. In Fig. 3b the length effect on the graphene thermal conductivity on the basis of MTP is compared with other full DFT studies, which show close trends. We note that depending on the type of exchange-correlation functional, supercell size and cut-off distance the estimated thermal conductivity of graphene by the change, which explains the remarkable scattering in the available data in the literature. On the basis of PBE/GGA functional and using the ShengBTE package, the thermal conductivity of monolayer graphene at room temperature was estimated to be 193647, 310048, 355049, 384550, 424251 and 5500 W/mK52. It is noticeable that MTP based estimate for the graphene lattice thermal conductivity agrees well with experimental and numerical results using the full DFT simulations. In Fig. 3c the phonon dispersion relation of graphene predicted by the DFT and MTP are compared, which reveal a remarkably close agreement. As it is well known, acoustic phonons are the main heat carriers in the graphene. In Fig. 3d we thus compare the contribution of different acoustic modes on the overall thermal conductivity by MTP and previous theoretical studies, which also show remarkable accuracy. Our results reveal that MTP potentials can be effectively used to estimate the lattice thermal conductivity not only via the NEMD simulations but also with the full iterative solution of the Boltzmann transport equation. Fig. 3, (a) NEMD estimations for the length effect on the room temperature lattice thermal conductivity of single-layer borophene along the armchair and zigzag directions (continuous lines illustrate the fits to the NEMD data points). (b) Lattice thermal conductivity of graphene at the room temperature as a function of mean free path by fully iterative solutions of the Boltzmann transport equation using the MTP (present study) and full DFT solutions by Fugallo et al. 53, Peng et al. 51, Gao et al. 54 and Qin et al. 47 with different exchange correlation functions. (c) Phonon dispersion relation of graphene acquired by the DFT (dotted line) and MTP (continuous line). (d) Contribution of ZA, TA and LA acoustic modes on the total lattice thermal conductivity of graphene by MTP (present study) and previous studies by Lindsay et al. 55, Gao et al. 54 and Qin et al.56. At this stage, we are now capable to explore the effective lattice thermal conductivity of graphene/borophene heterostructures by employing the FEM simulations, in which we employed ABAQUS/Standard along with python scripting. For the construction of heterostructures, we developed polycrystalline samples made of 5000 individual grains on the basis of Voronoi cells with mirror symmetry at all edges57. Then according to the volume 7 fraction, different grains were randomly assigned to be either graphene or borophene, simply by defining the corresponding thermal conductivity values acquired in the previous section. Since the borophene exhibits anisotropic thermal conductivity, for the corresponding cells the anisotropic thermal conductivity tensors were defined by randomly defining the orientation. The NEMD results for the thermal conductance of graphene/borophene grain boundaries were randomly chosen to introduce the interfacial conductance of every line connecting dissimilar crystals, and assuming perfect bonding (infinite conductance) for the rest of interfaces. To systematically investigate the size effect, we assume the equivalent grain size of the original polycrystalline sample as the domain size, assuming a square geometry for the equivalent average grain size57. A sample of heterostructure with 60% and 40% content of graphene and borophene phases, respectively, is shown in Fig. 4a. For the loading condition, we attached two highly conductive strips to the constructed sample and inward and outward the same value heat-fluxes (hf) were applied on the outer surfaces of these strips. As the initial value for the problem, the temperature of the outer surface of the cold strip (with outward flux) was set to zero. By solving the steady-state heat transfer problem, as shown in Fig. 4b a temperature profile establishes along the loading direction, which can be used to evaluate the effective thermal conductivity57. As expected and observable from the results shown in Fig. 4b-d, due to the high contrast in thermal conductivity of graphene and borophene, the temperature and heat flux profiles show highly non-uniform distributions. Such that the majority of heat fluxes are carried out via the percolation networks made of graphene crystals. To provide a more comprehensive vision on the heat transfer mechanism, we examined the effective lattice thermal conductivity for five different heterostructure compositions and for domain sizes ranging from 1 nm to 100 μm and the acquired results are illustrated in Fig. 4e-h. As the first impression, it is noticeable that for all samples with extremely large domain sizes around 100 μm, the effective thermal conductivity is not yet fully converged, revealing the importance of assuming the interfacial thermal conductance for the modelling of thermal transport in these heterostructures. The presented results reveal three main steps within the effective thermal conductivity with respect to the domain size. In the first step, occurring usually for domain sizes below 10 nm, the lattice thermal conductivity stays almost insensitive with respect to the domain size. This observation reveals that due to the presence of interfacial resistances, the embedded phases basically do not contribute to the heat flux transfer and they exhibit a void like behaviour. This issue is noticeable when comparing the thermal conductivity of heterostructures with 10% and 20% content of graphene nanosheets (find Fig. 4h), in which the sample with the higher content of the ultrahigh conductive crystals yields lower conductivity for domain sizes lower than 100 nm. The second step within the thermal transport mostly occurs for domain sizes from 10 nm to 10 μm, in which a smooth increase in the thermal conductivity by increasing the domain size is observable. Such a trend implies that the effect of interfacial resistance starts to decline by increasing the domain size. Let’s consider the sample with 10% content of borophene with a relatively large domain size of 250 nm, as observable from the results shown in Fig. 4c the borophene crystals contribute marginally to the heat flux transfer. For the sample with 40% 8 content of borophene, it is noticeable that the majority of heat flux is transferred by graphene networks percolating each other (find Fig. 4d). In this case, by increasing the content of borophene crystals they not only participate marginally in the heat transfer but also impede thermal transport within the highly conductive graphene grains. In the third and last step, which dominates the thermal transport for domain sizes larger than 10 μm, the thermal conductivity reaches a plateau and only slightly increases by a further increase of the domain size, which reveals that the effect of interfacial resistance starts to vanish. From a practical point of view, this second step of heat transfer would be more close to real experimental samples, because it is normally very difficult to make heterostructures with domain sizes larger than 0.01 cm. Our results for domain sizes from 10 nm to 10 μm highlight that within these ranges the interfacial thermal resistance plays the critical role and therefore should be taken into consideration. We would like to also clearly remind that in this study we mainly studied the lattice thermal conductivity, which may not be exactly the same as the total thermal conductivity, in which electrons contributions to the thermal conductivity should be also taken into account. Fig. 4, (a) A samples of constructed continuum model of graphene/borophene heterostructure with 40% content of borophene crystals to evaluate the effective lattice thermal conductivity of polycrystalline graphene structures, (b) established steady-state temperature profile for the same sample with domain sizes of 10 nm. (c and d) Samples of heat flux distributions. (e to h) Normalized effective thermal conductivity of heterostructure with respect to the graphene’s thermal conductivity. Conclusion In conclusion, our study confirms that machine-learning interatomic potentials trained over short ab-initio molecular dynamics trajectories enable efficient first-principles multiscale modeling via hierarchical employment of density functional theory/classical molecular dynamics/finite element simulations. In other words, it is possible to examine the properties and responses of novel complex microstructures, without prior knowledge concerning the 9 properties of building blocks. To show this novel possibility, we explored the lattice thermal conductivity of graphene/borophene heterostructures, a system that to the best of our knowledge there exists no viable classical modeling alternative. Furthermore, it is shown that the developed machine-learning interatomic potentials can be effectively employed to acquire the lattice thermal conductivity not only by classical molecular dynamics simulations but also with the full iterative solution of the Boltzmann transport equation. First-principles multiscale modeling is believed to offer novel and computationally efficient possibilities to evaluate the properties and improve the design of advanced nanostructures. Methods First-principles DFT calculations in this work were carried out using the Vienna Ab-initio Simulation Package (VASP)58–60. The generalized gradient approximation (GGA) and Perdew−Burke−Ernzerhof (PBE)61 functional was adopted in the calculations. We assumed the plane-wave cutoff energies of 500 eV in our simulations. The phonon dispersion and second- order force constants of graphene were obtained by density functional perturbation theory (DFPT) simulations over a 10×10×1 supercell sample using a 3×3×1 Monkhorst-Pack62 k-point grid along with the PHONOPY code63. Ab-initio molecular dynamics (AIMD) simulations were performed with a time step of 1 fs uing a 3×3×1 k-point gird. Acknowledgment B.M. and X.Z. 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). E.V.P and A.V.S. 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. References 1. Mounet, N. et al. Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds. Nat. Nanotechnol. 13, 246–252 (2018). 2. Oganov, A. R. & Glass, C. W. Crystal structure prediction using ab initio evolutionary techniques: principles and applications. J. Chem. Phys. 124, 244704 (2006). 3. Oganov, A. R., Lyakhov, A. O. & Valle, M. How evolutionary crystal structure prediction works-and why. Acc. Chem. Res. 44, 227–237 (2011). 4. Mannix, A. J. et al. Synthesis of borophenes: Anisotropic, two-dimensional boron polymorphs. Science (80-. ). 350, 1513–1516 (2015). 5. Feng, B. et al. Experimental Realization of Two-Dimensional Boron Sheets. Nat. Chem. 8, 563–568 (2016). 6. Zhou, X. F. et al. Semimetallic two-dimensional boron allotrope with massless Dirac fermions. Phys. Rev. Lett. 112, (2014). 7. Zhang, Z., Yang, Y., Gao, G. & Yakobson, B. I. Two-Dimensional Boron Monolayers Mediated by Metal Substrates. Angew. Chemie 127, 13214–13218 (2015). 8. Sun, B. & Barnard, A. S. Visualising multi-dimensional structure/property relationships 10 with machine learning. J. Phys. Mater. 2, 34003 (2019). 9. Oda, H., Kiyohara, S. & Mizoguchi, T. Machine learning for structure determination and investigating the structure-property relationships of interfaces. J. Phys. Mater. 2, 34005 (2019). 10. Schleder, G. R., Padilha, A. C. M., Acosta, C. M., Costa, M. & Fazzio, A. From DFT to machine learning: recent approaches to materials science–a review. J. Phys. Mater. 2, 32001 (2019). 11. Schmidt, J., Marques, M. R. G., Botti, S. & Marques, M. A. L. Recent advances and applications of machine learning in solid-state materials science. npj Comput. Mater. 5, 83 (2019). 12. Han, B. et al. Deep-Learning-Enabled Fast Optical Identification and Characterization of Two-Dimensional Materials. (2019). 13. Gubaev, K., Podryabinkin, E. V., Hart, G. L. W. & Shapeev, A. V. Accelerating high- throughput searches for new alloys with active learning of interatomic potentials. Comput. Mater. Sci. 156, 148–156 (2019). 14. Podryabinkin, E. V., Tikhonov, E. V., Shapeev, A. V. & Oganov, A. R. Accelerating crystal structure prediction by machine-learning interatomic potentials with active learning. Phys. Rev. B 99, 064114 (2019). 15. Ladygin, V. V, Korotaev, P. Y., Yanilkin, A. V & Shapeev, A. V. Lattice dynamics simulation using machine learning interatomic potentials. Comput. Mater. Sci. 172, 109333 (2020). 16. Novikov, I. S. & Shapeev, A. V. Improving accuracy of interatomic potentials: more physics or more data? A case study of silica. Mater. Today Commun. 18, 74–80 (2019). 17. Korotaev, P., Novoselov, I., Yanilkin, A. & Shapeev, A. Accessing thermal conductivity of complex compounds by machine learning interatomic potentials. Phys. Rev. B 100, 144308 (2019). 18. Podryabinkin, E. V & Shapeev, A. V. Active learning of linearly parametrized interatomic potentials. Comput. Mater. Sci. 140, 171–180 (2017). 19. Liu, X. & Hersam, M. C. Borophene-graphene heterostructures. Sci. Adv. 5, eaax6444 (2019). 20. Ghosh, S. et al. Extremely high thermal conductivity of graphene: Prospects for thermal management applications in nanoelectronic circuits. Appl. Phys. Lett. 92, 1–4 (2008). 21. Balandin, A. A. et al. Superior thermal conductivity of single-layer graphene. Nano Lett. 8, 902–907 (2008). 22. Jauregui, L. A. et al. Thermal Transport in Graphene Nanostructures: Experiments and Simulations. in 28, 73–83 (2010). 23. Cai, W. et al. Thermal transport in suspended and supported monolayer graphene grown by chemical vapor deposition. Nano Lett. 10, 1645–1651 (2010). 24. Tersoff, J. Modeling solid-state chemistry: Interatomic potentials for multicomponent systems. Phys. Rev. B 39, 5566–5568 (1989). 25. Stuart, S. J., Tutein, A. B. & Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 112, 6472–6486 (2000). 26. Brenner, D. W. et al. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J. Phys. Condens. Matter (2002). doi:10.1088/0953-8984/14/4/312 27. Lindsay, L. & Broido, D. A. Optimized Tersoff and Brenner empirical potential 11 parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene. Phys. Rev. B - Condens. Matter Mater. Phys. 81, 205441 (2010). 28. Wei, Z., Ni, Z., Bi, K., Chen, M. & Chen, Y. In-plane lattice thermal conductivities of multilayer graphene films. Carbon N. Y. 49, 2653–2658 (2011). 29. Hong, Y., Ju, M. G., Zhang, J. & Zeng, X. C. Phonon thermal transport in a graphene/MoSe2 van der Waals heterobilayer. Phys. Chem. Chem. Phys. 20, 2637– 2645 (2018). 30. Thomas, J. A., Iutzi, R. M. & McGaughey, A. J. H. Thermal conductivity and phonon transport in empty and water-filled carbon nanotubes. Phys. Rev. B 81, 045413 (2010). 31. Mortazavi, B. & Rabczuk, T. Multiscale modeling of heat conduction in graphene laminates. Carbon N. Y. 85, 1–7 (2015). 32. Fan, Z. et al. Thermal conductivity decomposition in two-dimensional materials: Application to graphene. Phys. Rev. B 95, (2017). 33. Hatam-Lee, S. M., Rajabpour, A. & Volz, S. Thermal conductivity of graphene polymorphs and compounds: From C3N to graphdiyne lattices. Carbon N. Y. (2020). doi:https://doi.org/10.1016/j.carbon.2020.02.007 34. Raeisi, M., Ahmadi, S. & Rajabpour, A. Modulated thermal conductivity of 2D hexagonal boron arsenide: a strain engineering study. Nanoscale 11, 21799–21810 (2019). 35. Bazrafshan, S. & Rajabpour, A. Engineering of thermal transport in graphene using grain size, strain, nitrogen and boron doping; a multiscale modeling. Int. J. Heat Mass Transf. 123, 534–543 (2018). 36. Shapeev, A. V. Moment tensor potentials: A class of systematically improvable interatomic potentials. Multiscale Model. Simul. 14, 1153–1173 (2016). 37. Liu, Y. & Yakobson, B. I. Cones, pringles, and grain boundary landscapes in graphene topology. Nano Lett. (2010). doi:10.1021/nl100988r 38. van der Zande, A. M. et al. Grains and grain boundaries in highly crystalline monolayer molybdenum disulphide. Nat. Mater. 12, 554–61 (2013). 39. Najmaei, S. et al. Vapour phase growth and grain boundary structure of molybdenum disulphide atomic layers. Nat. Mater. 12, 754–9 (2013). 40. Zhou, W. et al. Intrinsic structural defects in monolayer molybdenum disulfide. Nano Lett. 13, 2615–2622 (2013). 41. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 117, 1–19 (1995). 42. Mortazavi, B. et al. Strong thermal transport along polycrystalline transition metal dichalcogenides revealed by multiscale modeling for MoS 2. Appl. Mater. Today 7, 67–76 (2017). 43. Barrios-Vargas, J. E. et al. Electrical and Thermal Transport in Coplanar Polycrystalline Graphene–hBN Heterostructures. Nano Lett. 17, 1660–1664 (2017). 44. Schelling, P. K., Phillpot, S. R. & Keblinski, P. Comparison of atomic-level simulation methods for computing thermal conductivity. Phys. Rev. B 65, 1–12 (2002). 45. Zhang, X. et al. Thermal conductivity of silicene calculated using an optimized Stillinger-Weber potential. Phys. Rev. B - Condens. Matter Mater. Phys. 89, (2014). 46. Li, W., Carrete, J., Katcho, N. A. & Mingo, N. ShengBTE: A solver of the Boltzmann transport equation for phonons. Comput. Phys. Commun. 185, 1747–1758 (2014). 47. Qin, G., Qin, Z., Wang, H. & Hu, M. On the diversity in the thermal transport 12 properties of graphene: A first-principles-benchmark study testing different exchange-correlation functionals. Comput. Mater. Sci. (2018). doi:10.1016/j.commatsci.2018.05.007 48. Wang, H. et al. 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–2483 (2018). 49. Gao, Y. et al. First-principles study of intrinsic phononic thermal transport in monolayer C3N. Phys. E Low-Dimensional Syst. Nanostructures 99, 194–201 (2018). 50. Tan, X. et al. High thermoelectric performance in two-dimensional graphyne sheets predicted by first-principles calculations. Phys. Chem. Chem. Phys. 17, 22872–22881 (2015). 51. Peng, B. et al. Phonon transport properties of two-dimensional group-IV materials from ab initio calculations. Phys. Rev. B 94, (2016). 52. Taheri, A., Da Silva, C. & Amon, C. H. First-principles phonon thermal transport in graphene: Effects of exchange-correlation and type of pseudopotential. J. Appl. Phys. 123, 215105 (2018). 53. Fugallo, G. et al. Thermal conductivity of graphene and graphite: Collective excitations and mean free paths. Nano Lett. 14, 6109–6114 (2014). 54. Gao, Y. et al. First-principles study of intrinsic phononic thermal transport in monolayer C3N. Phys. E Low-dimensional Syst. Nanostructures 99, 194–201 (2018). 55. Lindsay, L., Broido, D. A. & Mingo, N. Flexural phonons and thermal transport in multilayer graphene and graphite. Phys. Rev. B 83, 235428 (2011). 56. Qin, Z., Qin, G., Zuo, X., Xiong, Z. & Hu, M. Orbitally driven low thermal conductivity of monolayer gallium nitride (GaN) with planar honeycomb structure: a comparative study. Nanoscale 9, 4295–4309 (2017). 57. Mortazavi, B., Pötschke, M. & Cuniberti, G. Multiscale modeling of thermal conductivity of polycrystalline graphene sheets. Nanoscale 6, 3344–3352 (2014). 58. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15–50 (1996). 59. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186 (1996). 60. Kresse, G. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 59, 1758–1775 (1999). 61. Perdew, J., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868 (1996). 62. Monkhorst, H. & Pack, J. Special points for Brillouin zone integrations. Phys. Rev. B 13, 5188–5192 (1976). 63. Togo, A. & Tanaka, I. First principles phonon calculations in materials science. Scr. Mater. 108, 1–5 (2015). 13