Unified deep learning network for enhanced accuracy in predicting thermal conductivity of bilayer graphene, hexagonal boron nitride, and their heterostructures

The thermal transport properties of layered materials have garnered considerable interest due to their pronounced anisotropy, which stems from the robust covalent bonds within the layers and the comparatively weak van der Waals interactions between them^1–5 This distinctive structure bestows upon these materials exceptional in-plane thermal conductivity coupled with a relatively diminished interlayer thermal conductivity, rendering them highly suitable for applications in electronic device heat dissipation. Graphene and hexagonal boron nitride (h-BN) are prominent among such materials, which excel in efficient in-plane thermal management and safeguard the underlying circuitry against the perils of overheating^6–10 The ease with which these materials can be exfoliated and subsequently stacked in rotational arrangements further gives rise to moiré superlattices, possessing captivating electronic and thermal attributes^5,11–22

The advent of moiré superlattices has precipitated an intensive wave of research focused on unraveling their unique electronic and thermal behaviors. Notable findings include the observation of Van Hove singularities in twisted graphene configurations, as well as the intriguing superconductor-to-insulator transition brought about by interlayer resonance enhancement.²³ Furthermore, the substantial reduction observed in thermal conductivity, alongside the marked enhancement in thermoelectric performance within moiré superlattices, has been associated with novel umklapp scattering processes tied to alterations in the Brillouin zone^19,24–26 Concurrently, van der Waals (vdW) heterostructures have emerged as a focal point of investigation, often exhibiting performances that either surpass or fall short of those exhibited by their constituent parent materials.^27–32 The capacity to modulate physical properties through straightforward stacking techniques is particularly alluring, and studies have revealed that the creation of vdW heterostructures significantly alters interlayer interactions, thereby profoundly influencing the thermal transport characteristics of the system. Examples include the temperature-dependent decrease in thermal resistance at the interface between graphene and two-dimensional SiC,³³ as well as the complex issue of thermal current convergence within graphene and h-BN heterostructures.^34,35

Additionally, thermal transport studies typically rely on precise computational frameworks such as the Density Functional Theory (DFT)³⁶and empirical potential-driven Molecular Dynamics (MD) simulations.³⁷ However, despite the high precision of DFT, its computational demands restrict the scale of systems that can be explored. On the other hand, while MD simulations accommodate larger systems, they often lack the comprehensive accuracy and applicability of DFT.^38,39 Moreover, investigations into moiré superlattice systems necessitate precise descriptions of the variations in interlayer coupling caused by twisting,^33,40–44 as well as addressing the increased atom count within the primitive cell and the concomitant escalation in computational costs due to reduced symmetry.^45,46 It is evident that a method capable of delivering DFT-level precision, while achieving the computational efficiency of MD simulations is required. Relying solely on DFT would result in prohibitively high computational expenses, whereas exclusive dependence on MD might not guarantee the necessary level of accuracy. Particularly when describing interlayer coupling, the conventional Lennard-Jones (L-J) potentials may prove inadequate to ensure veracity.^47,48 Thus, there is an urgent need for a novel approach that strikes a balance between computational precision and cost.

Recent advances in machine learning represent a promising solution, capable of achieving the accuracy of DFT, while retaining the computational efficiency of MD. Notable methodologies such as the Behler–Parrinello neural network (BPNN),⁴⁹ Gaussian Approximation Potential (GAP),⁵⁰ Gradient Domain Machine Learning Potential (GDMLP),⁵¹ Matrix Tensor Potential (MTP),⁵² Deep Learning Molecular Dynamics (DeepMD),⁵³ and Neural Evolutionary Potential (NEP)⁵⁴ represent key frameworks in this domain. For instance, the universal C element machine learning model, developed within the Gaussian Approximation Potential (GAP) framework,⁵⁵ along with the MLP crafted under the Neural Evolutionary Potential (NEP) framework specifically for layered h-BN,⁵ have both achieved near-DFT level accuracy, while substantially curtailing computational expenditures. Although MLP has been successfully applied to layered systems such as graphene and h-BN, its application within twisted graphene and h-BN heterostructures remains somewhat limited. This limitation stems from the intricate potential energy landscapes of heterostructures’ twisted systems and the feeble long-range van der Waals forces, which collectively pose significant challenges in devising appropriate machine learning potential functions.^56,57

In this research, we employed DFT to conduct rigorous ab initio molecular dynamics (AIMD) simulations and detailed static calculations on graphene, h-BN, and their complex heterostructures. These materials were carefully subjected to a range of mechanical strains, perturbations, twist angles, and defects to gain a comprehensive understanding of their behavior. Utilizing the wealth of energy, force, and virial data obtained from these simulations, we constructed a meticulously curated training set encompassing 1253 distinct structural configurations. By harnessing the power of the NEP framework, seamlessly integrated within the Graphics Processing Units Molecular Dynamics (GPUMD) environment, we developed a highly accurate machine learning potential (MLP) model. This model faithfully mirrored the intricate potential energy surface predicted by DFT, ensuring both precision and reliability in our subsequent analyses. Thorough validation procedures, including the meticulous examination of lattice constants and phonon dispersion relations, further confirmed the robustness of our MLP model. This solid foundation enabled us to further explore the intricate thermal transport properties of these two-dimensional materials and their heterostructures. Our in-depth analysis revealed a profound impact of twist angles on thermal conductivity, particularly evident in heterostructures where a reduction exceeding 35% was observed. This substantial decline can be attributed to the twist angle-induced softening of phonon modes, coupled with a concurrent surge in phonon scattering rates. These phenomena collectively intensify the anharmonic interactions among phonons, significantly altering their transport properties. Our research not only underscores the remarkable efficacy of MLP models in elucidating the intricate thermal transport behavior of two-dimensional materials and their heterostructures but also provides valuable insights into the fundamental micro-mechanisms governing the influence of twist angles on thermal conductivity. In contrast to previous studies such as the universal Si MLP developed by Li et al.,⁵⁸ our research methodology is meticulously tailored to accommodate the unique properties of CBN elements and subsequent structures like twisted graphene and hexagonal boron nitride heterostructures. By adopting the Neural Evolution Potential (NEP) framework for constructing a unified descriptor, our approach aims to capture the complex interactions present in these systems, thereby broadening the scope of problems that our MLP can effectively address. This innovative approach offers a promising avenue for the rational design of next-generation thermal management materials.

In our study, we employed first-principles computations using the Vienna Ab Initio Simulation Package (VASP),⁵⁹ a powerful tool based on DFT.⁶⁰ This allowed us to accurately model the complex dynamics between ions and electrons, which are essential to our system. We achieved this using the Projector Augmented Wave (PAW)⁶¹ method, in conjunction with an advanced plane wave basis set.

To faithfully depict the intricacies of the exchange-correlation energy terms, we integrated the Perdew–Burke–Ernzerhof (PBE) function^62,63within the framework of the Generalized Gradient Approximation (GGA).^64,65 In order to correctly identify the van der Waals interaction, the vdW-DFT-cx method (optB86-vdW)⁶⁶ was used in the calculation. The Brillouin zone was sampled using a Γ-centered grid with a linear KSPACING spacing of approximately 0.25 Å⁻¹ and Gaussian smearing with a width of 0.01 eV. For force calculations, a finer support grid was utilized to enhance numerical accuracy. All calculations were performed with a plane-wave energy cutoff of 520 eV. All parameters were established following the previous study.⁵ 

The foundation of our investigation lies in the construction of a comprehensive and diverse dataset comprising 1253 structures for graphene/h-BN systems (see details in Table S1 in the supplementary material). Delving into the intricacies of stacking arrangements and linking modes, our dataset explored a plethora of possibilities for graphene, h-BN, and their heterostructures. Recognizing the impact of larger unit cells in capturing nuanced structural nuances, we systematically considered the influence of external factors such as strain, defects, and perturbations. The resultant training set, hosting a staggering 165704 atoms, stands as a testament to our commitment to comprehensively explore the rich landscape of structural configurations.

We are dedicated to developing an accurate graphene/h-BN potential with the aid of utilizing the fourth-generation NEP from the GPUMD software package.^54,67 It entailed the strategic use of a multi-layer perceptron architecture with a meticulously designed single hidden layer. The intricacies of the atomic environment descriptor were brought to life by a judicious combination of Chebyshev basis functions for the radial part and Legendre polynomials for the three-body angle part. A carefully chosen cutoff of 8 Å for the radial component and 3.5 Å for the angular component ensured the fidelity of our MLP architecture, which consistently featured 50 neurons in the hidden layer throughout the entire training process.

Our pursuit of a comprehensive understanding extended beyond the confines of the theoretical framework, integrating machine learning potentials with various computational tools. Collaborating seamlessly with GPUMD, phonopy, and shengBTE,^54,68,69 we embarked on an exploration to extract phononic properties and correlate them with the intricate thermal behaviors of our system.

The process of constructing a complete MLP and predicting the thermal transport properties through MLP is shown in Fig. 1. We search for training structures to obtain energy, force, and virial information for machine learning to acquire MLP, and get the corresponding force constants through the interface with different software and MLP. We then obtain thermal transport-related information using shengBTE, a software that solves the phonon Boltzmann transport equation. In the construction of the training set, we start with pristine multilayer graphene and h-BN as the initial structures, simulating 5000 steps with an AIMD time step of 2 fs, taking one configuration every 100 steps for a constant pressure annealing simulation at 300–1500 K. After obtaining the initial configurations, we apply perturbations, strains, and defects to the initial configurations and perform static calculations to obtain the corresponding energy, force, and virial information to better fit an appropriate potential energy surface. The 2–4 structures were randomly chosen from the displacement POSCAR files generated using phonopy. These selected configurations served as the perturbed configurations, augmenting the diversity within our dataset and strengthening the robustness of our analysis. Moreover, larger supercells have more positional and angular information, providing more operational space for building an MLP suitable for twisted systems. Larger supercells offer more operational space for constructing MLP for distorted systems like twisted graphene, accommodating varied atomic arrangements and capturing intricate structural features during training. This enhancement in structural detail is crucial for accurately modeling bonding and interaction phenomena in distorted systems, enriching the robustness and fidelity of the study.^70,71 Essentially, twisted systems are caused by the deviation of the original atomic positions and angles, which changes the potential energy surface. If all possible positions and angles can be expressed, information on different twist angles will also be included. Previous research found that a larger rotating crystal cell system contains rich positional and angular information. Therefore, we constructed twisted structures of graphene and h-BN, with no fewer than 100 atoms in each possible twisted structure. We further choose structures with twist angles of 21.78°, 27.8°, 13.17°, and 9.43° from among these structures for AIMD annealing simulations at 300–1500 K to observe the evolution of twist angles at different temperatures. In our AIMD simulations, conducted over the temperature range of 300 to 1500 K, we replicated an annealing process to emulate real-world conditions. Throughout these simulations, we saved configurations every 100 fs, allowing for a comprehensive sampling of thermal fluctuations. This meticulous approach ensures the robustness of our training dataset.

Notably, our trained MLP meticulously encompasses essential features of multilayer systems, including a range of stackings, interlayer distances, and structural arrangements. By leveraging DFT data calculated with VASP as the foundation of our training set, we introduce an intrinsic periodicity, reinforcing the MLP’s capacity to accurately articulate the properties of multilayer systems^72,73 (for further details, refer to Fig. S2 in the supplementary material).

Constructing a unified potential for graphene, h-BN, and their heterostructures is very challenging. On the one hand, it needs to consider all possible situations of graphene and h-BN alone, which has been solved in the initial construction of the training set; on the other hand, the problem of lattice mismatch in heterostructures becomes the key to restricting the unified potential. However, graphene and h-BN have very similar lattice constants, making their layered heterostructures naturally have low lattice mismatch, which solves the second problem smoothly. Therefore, we constructed two initial structures of graphene/h-BN lateral heterostructures and vdW heterostructures and considered twist angles in vdW heterostructures as in single graphene and h-BN to fully consider possible atomic positions and angles. Finally, we constructed a training set containing 1253 structures with 165704 atoms, fully obtaining energy, force, and virial information to realize the scanning of potential energy surfaces in machine learning processes and performing MLP of graphene, h-BN, and their heterostructures under a unified descriptor.

Through the fourth-generation NEP framework built into GPUMD, we fit the training set with MLP to obtain a unified MLP. In this process, we set the training steps to 350 000 steps to obtain a convergent MLP. As shown in Fig. 2(a), when the training steps reach 200 000 steps, the loss functions of energy, force, and virial have fully converged, indicating that a fully convergent MLP has been obtained. A further comparison of the trained energy, force, and virial with DFT results is shown in Figs. 2(b)–2(d). The root mean square errors (RMSEs) for energy, force, and virial are 1.8 meV/atom, 78 meV/Å, and 51 meV/atom, respectively. Comparatively, previously reported RMSE values for h-BN are 2 meV/atom, 45 meV/Å, and 13 meV/atom, while for graphene, they are 1 meV/atom, 43 meV/Å, and 8 meV/atom.^{5, 55} These results indicate that the trained MLP exhibits high consistency with DFT and achieves DFT-level accuracy.

While maintaining a high consistency in energy, force, and virial with DFT during training tests is important, it does not guarantee similar consistency during the application process, especially in the study of thermal transport. In this context, an accurate prediction of lattice constants and dispersion relationships becomes essential. It is widely acknowledged that the precision of lattice constants and dispersion relationships significantly influences the ability to predict lattice thermal conductivity within a reasonable order of magnitude. Therefore, we verified the lattice constants of graphite, graphene, AB stacked graphene, bulk h-BN, monolayer h-BN, and AB stacked h-BN compared with experiments or DFT calculations. As shown in Table I, where “a” and “c” represent the lattice directions, CBN-NEP shows good predictive performance in all verified systems. This suggests that CBN-NEP can correctly identify lattice constants.

In the phonon gas model, the depiction of lattice vibration relies on the spacing between atoms. CBN-NEP’s precise prediction of lattice constants establishes a strong basis for accurately forecasting lattice vibration. The accuracy of identifying lattice vibration can be assessed through the phonon dispersion relationship, by comparing the results of DFT and CBN-NEP calculations for AB stacked graphene and h-BN, as depicted in Figs. 3(a) and 3(b). It becomes evident that CBN-NEP excels in capturing information related to lattice vibrations, particularly in the phonon dispersion relationship. Unlike traditional empirical potential fields that lack accuracy in describing high-frequency phonons, CBN-NEP demonstrates excellent performance in characterizing high-frequency lattice vibration.

In addition to harmonic lattice vibration information in phonon transport, anharmonic interactions among phonons are also crucial. The primary approach for applying DFT to phonon transport research involves obtaining the harmonic and anharmonic interaction force constants between atoms, known as the second-order (IFC2) and third-order (IFC3) force constants. We use the second-order and third-order force constants as input parameters to solve the phonon Boltzmann transport equation, extracting phonon transport-related information. This involves integrating CBN-NEP with phonopy and thriorder.py to compute the force constants of AB stacked graphene. Then, using the shengBTE software package, we derive the corresponding thermal conductivity. In Fig. 3(c), we show the thermal conductivity acquired using both DFT and CBN-NEP from 300 K to 1500 K. A comparison of thermal conductivity demonstrates the consistency of CBN-NEP with DFT across a broad temperature spectrum. The thermal conductivity at room temperature for AB stacked graphene and h-BN is 1893 and 246 W/m-K, respectively, aligning with those of previous research findings.^10,79–84 Additionally, the Grüneisen parameter serves as an indicator of the strength of anharmonic interactions in nonharmonic interactions, making it a pivotal parameter for nonharmonicity assessment. Therefore, we compare the total Grüneisen parameter between DFT and CBN-NEP results, as illustrated in Fig. 3(c). Across the temperature range of 300 to 1500 K, both DFT and CBN-NEP obtained total Grüneisen parameters exhibit consistent trends. This indicates that CBN-NEP effectively discerns harmonic and nonharmonic vibration information.

Up to now, we have demonstrated the consistency of CBN-NEP with DFT in bilayer graphene and h-BN concerning phonon dispersion relations, lattice thermal conductivity, and anharmonicity. Next, we extended the application of CBN-NEP to explore the thermal behavior of double-layer twisted graphene, h-BN, and graphene/h-BN heterostructures, affirming its suitability for twisted systems. In these systems, the interlayer twist angle disrupts the original symmetry of the structure, leading to a dramatic increase in the number of atoms within the unit cell under twisting. For instance, previous research on twisted graphene has indicated that a twist angle of 21.78° results in the smallest unit cell beyond the AA and AB stacking configurations.^45,46 Limited by the memory requirements of solving the Boltzmann transport equation, we focus solely on the 21.78° twist structure for research, with all subsequent twist angles referring specifically to 21.78°unless otherwise stated.

Through the integration of CBN-NEP and shengBTE software packages, we examined the thermal transport properties of graphene, h-BN, and graphene/hexagonal boron nitride heterostructures at a specific twist angle. As depicted in Figs. 4(a)–4(c), the absence of imaginary frequencies for these materials under twisted configurations signifies their structural stability at the given angle. Notably, within the twisted systems, there is a noticeable softening of low-frequency phonons, with the out-of-plane acoustic (ZA) mode near the gamma point exhibiting an almost flatband profile. This phenomenon leads to intensified phonon scattering and coupling, corroborating findings from earlier research endeavors.

The previous overlap of optical phonons and acoustic phonons influences the behavior of phonon scattering, which is governed by the conservation of energy and momentum, thus impacting the material’s thermal conductivity. As illustrated in Fig. 4(d), the calculated thermal conductivity values for T graphene and T h-BN at 300 K are 1530 and 123 W/m-K, respectively. Furthermore, the thermal conductivity of T graphene/h-BN heterostructures is computed at 305.04 W/m-K. In comparison with the AB-stacked heterostructure’s thermal conductivity, the thermal conductivity of the twisted graphene/boron nitride heterostructures has plummeted by over 35%, amounting to merely one-fifth of the twisted graphene's value, yet still significantly surpassing that of the twisted h-BN. The phonon dispersion, thermal conductivity, and Grüneisen parameter of the AB stacked graphene/h-BN heterostructure, calculated by using both DFT and NEP, are illustrated in Fig. S1 in the supplementary material. The results from both approaches demonstrate consistency, underscoring the accuracy of the CBN-NEP potential. Furthermore, both graphene and h-BN exhibit a decline in their twisted thermal conductivities compared with their AB-stacked counterparts, aligning with previous observations on the diminishing effect of twisting on thermal conductivity.

Delving deeper into the underlying physical principles, we have analyzed the critical equation determining thermal conductivity: $k= 1 3Σ C V v 2τ$⁠, where $C V$ denotes the phonon specific heat capacity, v represents the phonon group velocity, and $τ$ stands for the phonon lifetime, inversely related to the phonon scattering rate. Figures 5(a)–5(c) illuminate the variations in phonon group velocity, indicating a substantial reduction due to twist-induced phonon mode softening and increased optical phonon–acoustic phonon overlap, which are particularly evident in the twisted heterostructures. Such deceleration often stems from altered interlayer interactions that significantly affect phonon–phonon interaction dynamics.

Additionally, we investigate the consequences of alterations in phonon scattering rates. As shown in Figs. 5(d)–5(e), transitioning from an AB-stacked to a twisted structure results in a marked elevation in the phonon scattering rate, indicating that twisting exacerbates anharmonic effects through modified interlayer interactions. The visible change in phonon scattering rates manifests this effect: enhanced anharmonicity leads to increased phonon scattering likelihood and curtailed phonon lifetimes. Within this context, the changes in the graphene system are relatively subtle in comparison with those in the h-BN and heterostructure systems, explaining the comparatively minor impact of twisting on graphene’s thermal conductivity. Conversely, the heterostructure system undergoes an order-of-magnitude shift in its scattering rate, precipitating a drastic diminution in phonon lifetimes and consequently impeding phonon transmission. Thus, the concurrent reduction in phonon group velocity and scattering rate culminates in a significant decrease in thermal conductivity under twisted conditions as opposed to the AB-stacked configuration. Despite the twist-induced increase in phonon mode overlap, homostructures do not exhibit the same extent of reduction in phonon group velocity and scattering rate as observed in heterostructures. These divergent outcomes ultimately lead to a more substantial reduction in thermal conductivity for heterostructures relative to homostructures when subjected to identical twist conditions.

We constructed a training set containing 1253 structures and 165704 atoms, which adequately captures energy, forces, and virial information. Through potential energy surface scanning in the machine learning process, we built a unified MLP for graphene, h-BN, and their heterostructures under a consistent descriptor, which we named CBN-NEP. The trained MLP shows high consistency with DFT and achieve DFT precision. Furthermore, we used CBN-NEP to investigate the thermal properties of bilayer twisted graphene, h-BN, and graphene/h-BN heterostructures to verify their applicability in twisted systems. The results indicate that the thermal conductivity of twisted graphene/boron nitride heterostructures drops by more than 35%, only one-fourth of that of twisted graphene, but still significantly exceeds the thermal conductivity of twisted h-BN. We conclude that the developed CBN-NEP is a promising tool for studying heat transport in materials with strong phonon anharmonicity or spatial disorder, which are often not accurately treated by traditional empirical potentials or perturbative methods. The MLP offers unparalleled accuracy and adaptability in detecting interlayer heat transfer phenomena in van der Waals materials, providing an innovative analytical tool for delving into the complexities of two-dimensional materials and their structures.

See the supplementary material for a comprehensive overview of the calculation details pertaining to the thermal conductivity, including the phonon dispersion of the three layers of graphene and h-BN, using both DFT and NEP methods. Additionally, the mapping of different elements to corresponding atomic numbers and structural counts is also contained.

This research was funded in part by the National Natural Science Foundation of China (Grant No. 12105242) and by the Yunnan Fundamental Research Project (Grant Nos. 202201AT070161, 202301AW070006, and 202401AT07031).

The authors have no conflicts to disclose.

R.C. and Y.T. contributed equally to this work.

Rongkun Chen: Conceptualization (equal); Data curation (equal); Investigation (equal); Writing – original draft (equal). Yu Tian: Formal analysis (equal); Methodology (equal); Validation (equal); Writing – original draft (equal). Jiayi Cao: Investigation (equal); Software (equal). Weina Ren: Conceptualization (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Shiqian Hu: Conceptualization (lead); Funding acquisition (lead); Supervision (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead). Chunhua Zeng: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.