First-principles-based machine learning interatomic potential for molecular dynamics simulations of 2D lateral MoS₂/WS₂ heterostructures

Two-dimensional (2D) transition-metal dichalcogenides (TMDs) exhibit excellent potential for applications in electronics, valleytronics, emerging quantum optoelectronics, and sensors^1–4 due to their unique electronic, mechanical, optical, and thermal properties. Monolayer TMDs, such as MoS₂, WS₂, and their heterostructures, have been successfully fabricated for nano-electromechanical system (NEMS) devices,^5,6 light-emitting diodes, high-electron mobility transistors (HEMTs),^7–11 etc. In these applications, materials are subject to more complex mechanical loads and Joule self-heating. As devices continue to decrease in size, the interfacial thermal resistance (ITR) of the heterostructures becomes a primary contributor to overall thermal resistance.^12,13 It is particularly relevant for lateral TMD–TMD′ heterostructures,^14–16 which are not only pivotal due to their contribution to thermal management but also serve as functional components in electronic devices.^17,18 Therefore, it is critical to investigate interfacial thermal conductance and mechanical properties of lateral TMD–TMD′ heterostructures to advance manufacturing processes and enhance the performance of these materials in diverse applications.

There are two conventional methods to present the atomic interactions in the structures and heterostructures, which is crucial for the study of their physical properties. Based on density functional theory (DFT),^19,20 ab initio molecular dynamics (AIMD)²¹ simulations can accurately describe structures and atomic trajectories, albeit at the expense of efficiency. Due to the consideration of complex electronic contributions, these simulations are typically limited to tens of picoseconds and hundreds of atoms. However, intricate scenarios, such as the heterostructures of TMDs, need longer time and larger cell sizes. To address these limitations, molecular dynamics (MD) simulations use the empirical potentials,^22–24 such as Stillinger−Weber (SW) potentials and reactive empirical bond-order (REBO) potentials, to exclusively describe atomic interactions. This enables the study of mechanical and thermal properties of structures on a larger scale (hundreds of nanoseconds and nanometers).^25,26 Nevertheless, the accuracy of empirical potential-based MD simulations heavily relies on the accuracy of empirical potentials. On the one hand, the empirical potentials usually overlook intricate details of atomic interactions, leading to an inability to precisely represent higher-order terms such as forces and virials, which are critical aspects for exploring the physical properties of micro–nano structures. On the other hand, the limited transferability of empirical potentials raises concerns, as different potentials may have varying scopes of application. Hence, it is usually questionable to use empirical potentials to calculate diverse properties such as ITR and Young's modulus of heterostructures simultaneously. Moreover, since certain properties are highly correlated with forces and virials, empirical potentials with fixed, simplistic forms struggle to accurately describe ITR of heterostructures and mechanical properties such as Young's modulus simultaneously. Therefore, there is an urgent need for a method capable of representing atomic trajectories accurately and efficiently.

As it stands, AIMD is inefficient when dealing with large-scale atom computations, while the accuracy of empirical potential-based MD simulations largely depends on the quality of potentials. A potential balancing the efficiency and precision is urgently required to ensure effective computational simulations. However, developing interatomic potentials involving fitting pre-selected function forms with ab initio energy surfaces has proven challenging due to inherent complexity.²⁷ Fortunately, with the development of artificial intelligence (AI) and the increasing application of machine learning (ML) in material science, the above challenges have been well addressed.^28,29 Various machine learning interatomic potentials (MLIPs) have been developed, including the Gaussian approximation potential (GAP),^30–32 TensorMol,³³ the neuroevolution machine learning potential (NEP),³⁴ and the deep neural network potential (NNP).³⁵ Notably, neural networks (NNs) possess the unique ability to theoretically represent unknown multidimensional real-valued functions with arbitrary precision by selecting appropriate network models. Specifically, NNPs can accurately reproduce the relationship between the nuclear positions and energies.³⁶ Recent research³⁷ demonstrates that the NNP model matches the accuracy of quantum mechanics for both finite and extended systems, showcasing its size-extensive nature. NNPs overcome the limitations associated with auxiliary quantities and introduce a flexible family of loss functions, exhibiting excellent performance. At present, MLIP has been successfully applied in diverse materials, including TiO₂,³⁸ CdTe,³⁹ Ga₂O₃,⁴⁰ and N–Ga–Al system,⁴¹ making the NNP with the broadest range of applications and the most promising prospects in current research.

In this work, an NNP for monolayer MoS₂/WS₂ and their heterostructures was developed, using the deep potential molecular dynamics (DeePMD)⁴² simulations and DP-Gen (Deep Generator),⁴³ achieving accuracy comparable to first-principles calculations. Subsequently, the obtained results were systematically compared with AIMD to validate the effectiveness of the MLIP we developed. Finally, MD using the MLIP was used to predict Young's modulus and thermal conductivity of the monolayer MoS₂/WS₂ and ITR of their lateral heterostructures. The paper is organized as follows: Sec. II provides a detailed account of the preparations of the initial dataset and the training of MLIP. Section III assesses the performance of MLIP through testing, and rigorous verification of its accuracy of MLIP is well verified. Subsequently, we present the calculated Young's modulus and thermal properties of the monolayer MoS₂/WS₂. Finally, a summary is given in Sec. IV.

To develop an MLIP, a dataset was initially curated that included snapshots of monolayer MoS₂/WS₂ and their heterostructures, as shown in Fig. 1. These structures were constructed by AIMD simulations. The lattice parameters of unit cell MoS₂ (a = 3.190 Å, b = 5.526 Å) and WS₂ (a = 3.184 Å, b = 5.515 Å) were established. The supercells, each consisting of N = 90 particles, were constructed for MoS₂ and WS₂. The supercell parameters of MoS₂ were set at a = 15.952 Å and b = 16.577 Å. Similarly, WS₂ was constructed in the same way. To enhance the precision of Young's modulus expression and facilitate MLIP training, stretch snapshots using SW potential⁴⁴ were also incorporated into the datasets. The materials underwent pre-stretched at a constant strain rate of 10⁹ to reach a 30% strain. The high strain is widely used in empirical potential-based MD simulations for calculating Young's modulus.^45,46 To build MoS₂/WS₂ heterostructures, W atoms were replaced by Mo atoms at a rate of 10%, resulting in structures such as Mo_0.1W_0.9S, Mo_0.2W_0.8S, Mo_0.3W_0.7S, Mo_0.7W_0.3S, Mo_0.8W_0.2S, and Mo_0.9W_0.1S. Specifically, the heterostructures with sharp interfaces along zigzag and armchair directions are included in MLIP's datasets to enhance the MLIP's capability to handle heterostructures. Random perturbations in the range of 0–0.05 Å were applied to the atomic positions of MoS₂/WS₂ and heterostructures to account for the effects of thermal vibrations, further enriching the dataset.

Atomic coordinates were acquired using Vienna Ab initio Simulation Package (VASP)⁴⁷ based on DFT to generate the initial datasets for MLIP. The projector-augmented wave (PAW) potentials were employed to describe interactions between electrons and ions, and the exchange-correlation functional of Perdew–Burke–Ernzerhof (PBE) was used to describe the electron exchange and correlation potential in all DFT calculations,^48–50 and it is widely used for the investigation of 2D TMDs including Young's modulus, thermal conductivity, etc., which agree well with experimental results.^51–57 Given the expensive scale of the simulation, only the Γ point in the Brillouin zone was sampled, with an energy cut-off set at 500 eV. The datasets contained atoms coordinate and the corresponding total energy of systems, interatomic forces, and virials of 34 000 datasets spanning the temperatures from 200 to 400 K. This compilation consists of 22 000 sets of MoS₂/WS₂ crystal structures within 200–400 K, 2000 sets of stretch structures at 300 K, and 10 000 sets of heterostructures with a sharp interface within 200–400 K. To ensure the robust performance of the model, the datasets were randomly divided into training and testing sets at a ratio of 5:1. This careful partitioning supports the reliability and generalizability of the NNP in accurately predicting the properties of the MoS₂/WS₂ structures across a diverse range of conditions.

The MLIP for monolayer MoS₂/WS₂ and their heterostructures was developed using the Deep Potential Molecular Dynamic package (DeePMD-kit)^{41,42,58–60} and DP-Gen (Deep Generator). To enhance the MLIP performance in intricate scenarios, DP-Gen was used to generate structures with fewer types in the dataset. A multilayer neural network was used to predict atom corresponding energies. One challenge often encountered in MLIP development is finding a suitable input configuration that preserves translational, rotational, and permutational symmetries of the environment.⁶¹ In this work, Smooth Edition of Deep Potential (DeepPot-SE)⁶² was used to transform structural descriptors, aligning with the aforementioned requirements. To obtain the structural information while ensuring the model's continuous differentiability, the cut-off radius r_c and smooth cut-off parameter r_cs were set at 7 and 6.50 Å, respectively. The activation function of NN is hyperbolic tangent. In addition, the hyper-parameters of the NN are shown in Table I.

A large-scale atomic/molecular massively parallel simulator (LAMMPS)⁶⁵ was used to validate MLIP's accuracy and investigate the mechanical and thermal properties of monolayer MoS₂/WS₂ and their lateral heterostructures. MLIP-based MD was carried out to mirror the simulation process conducted with AIMD, ensuring consistency between MLIP calculations and AIMD benchmarks. Subsequently, large-scale atomic simulations were executed to calculate thermal conductivity, Young's modulus, and ITR.

To validate the accuracy of MLIP, the results obtained from MLIP-based MD simulations were compared with those from AIMD simulations. The trained MLIP was employed to calculate energies, forces, and virials. Figure 2(a) shows the difference in dataset energies between MLIP-based MD and DFT-based AIMD simulations. Figure 2(b) illustrates the variation in atom forces along the x, y, and z directions. It is obvious that energies and forces consistently align along the diagonal between MLIP predictions and AIMD benchmark. The root mean square of error (RMSE) for energies and forces is remarkably low, at 3.48 meV per atom and 0.042 eV/Å, respectively, which is significantly lower than reported values in literature for MLIP-based MD.^{41,58,59,66,67} The maximum error of energies and forces between the MLIP prediction values and DFT-calculated results is 6.91 meV per atom and 0.24 eV/Å, respectively.

For large-scale molecular dynamics simulations, it is crucial to ensure the continuity and stability of the potential used. To assess this, changes in the potential energy and kinetic energy of the system were examined under a microcanonical (NVE) ensemble by MLIP-based MD simulations. A MoS₂/WS₂ lateral heterostructure was considered, containing 864 atoms with 16 × 9 unit cells with elements Mo₁₄₄S₅₇₆W₁₄₄, and it possessed supercell lattice parameters a = 51.05 Å and b = 49.73 Å. Initially, a periodic boundary condition was applied in the x direction, while a free boundary condition was applied in the y direction to represent monolayer MoS₂/WS₂ heterostructures. Subsequently, the NPT ensemble was used to equilibrate the system at 3 K and zero pressure for 50 ps, aiming to reduce the stress of the structures. Finally, the energy data were collected during the period of the NVE ensemble spanning 100 ps. Figure 3 indicates the trends of total energy, kinetic energy, and potential energy. It is evident that the total energy of the system, equal to the sum of potential energy and kinetic energy, remains stable. This observation provides unequivocal evidence that the MLIP developed is continuous and stable, affirming its suitability for conducting molecular dynamics simulations. In addition, the phonon dispersion was calculated using MLIP and compared with the results from experiments⁶⁸ and DFT calculations.^69,70 Here, the experimental results are from bulk MoS₂ since few experiments have been carried out to measure the phonon dispersion of monolayer MoS₂.⁶⁸ As shown in Fig. 4, the results of phonon dispersion using MLIP agree well with the DFT calculations.

The radial distribution functions (RDFs) g(r) describe the probability density concerning the distance from a reference particle. To assess the capability of MLIP in replicating atomic trajectories with accuracy comparable to DFT, we compared the RDFs calculated by DFT-based AIMD simulations and MLIP-based MD simulations. The chosen distance r for RDF is set to 8 Å, slightly larger than the cutoff r_c = 7.5 Å of MLIP. As shown in Fig. 5, the g(r) of different elements at 300 K demonstrates a remarkable alignment between the RDF predicted by MLIP and that calculated by DFT, indicating the excellent reliability of the constructed MLIP for the theoretical simulations of MoS₂/WS₂ and their lateral heterostructures.

In the fabrication of flexible devices, TMDs are often subjected to mechanical loads, emphasizing the need for a thorough analysis of their mechanical properties. As a representative example, Young's modulus was calculated to assess the accuracy of MLIP in predicting the mechanical properties of TMDs. Young's modulus is defined as the ratio of the stress applied to the resulting axial strain within the linear elastic region of the material, and it serves as a critical parameter in characterizing materials’ ability to resist deformation.

To verify the applicability of MLIP for the mechanical properties of TMDs, LAMMPS was used to simulate the stretching of TMDs. Various structures were built, incorporating MoS₂ and WS₂ along the armchair and zigzag direction, respectively. The supercell parameter is 50 × 50 Ų. The time step of 1 fs and Nose–Hoover thermostat were used in simulations. The NPT ensemble was employed to equilibrate the system at 300 K and zero pressure for 30 ps. Subsequently, the NPT ensemble was applied in the non-stretch direction at 300 K and zero pressure, and deformed materials in the tensile direction with a strain rate of 10⁹ to reach a 30% strain.

We calculated Young's modulus of MoS₂ and WS₂ as 223.5 and 246.5 GPa, respectively. Figure 6(a) illustrates the trend of the stress–strain curve in the zigzag direction of MoS₂ and WS₂. It reveals that MoS₂ possesses a lower Young's modulus compared to WS₂. The slope of the linear portion before a strain of 5% represents Young's modulus. Beyond a strain of 15%, the material structure becomes increasingly unstable, eventually leading to fracture at a strain of 25%. Figure 6(b) illustrates Young's modulus of the heterostructure when stress is applied along and across the interface, and Young's modulus of the heterostructure is 228 and 55 GPa, respectively. Notably, once stress is applied along the interface, the structure breakage occurs in the MoS₂ region at 25% strain, exhibiting Young's modulus of 228 GPa, resembling reminiscent of the behavior observed in MoS₂. The slightly reduced Young's modulus in the heterostructure may be attributed to stress concentration at the fracture site in the MoS₂ region, resulting in a highly unstable structure and eventual material fracture. Interestingly, a significantly smaller Young's modulus is observed, when stress is applied across the interface, most possibly due to unbalanced stress on sulfur atoms at the linkage of the interface, leading to the breakage of Mo–S bonds.

The MLIP-based calculations of Young's modulus in the armchair and zigzag directions for MoS₂ and WS₂ are shown in Table II, alongside available experimental data,^71–74 empirical potential-based MD simulations,⁷⁵ and DFT calculations.^55,76,77 The comparison highlights that the difference between the calculated value of MLIP-based MD and values reported in the literatures fall within a reasonable range, certifying that MLIP-based MD exhibits comparable accuracy to DFT.

A monolayer crystal structure was constructed for MoS₂ and WS₂ with the size of 60.6 × 60.8 Ų. Periodic boundary conditions are applied in the in-plane direction to mimic the infinite size of the structures. The Nose–Hoover thermostat was used in thermal conduction's calculation of NPT and NVT ensembles. The NPT ensemble was applied to equilibrate the system at 3 K and zero pressure for 200 ps with a time step of 1 fs. Then, the NVT ensemble was performed at 300 K. Following a 200 ps equilibration period in the NVT ensemble, the ensemble was switched to an NVE ensemble of 2 ns to collect data on the heat flux at 300 K.

Since the thermal conductivity of MoS₂ and WS₂ is isotropic in the in-plane direction, the values of thermal conductivity were average in the in-plane direction. Multiple simulations with different initial conditions were performed to reduce random errors, and the values of thermal conductivity were averaged to obtain the thermal conductivity κ of MoS₂ and WS₂ at 300 K. The thermal conductivity of MoS₂ and WS₂ is 105 ± 18 and 87 ± 15 W m⁻¹ K⁻¹, respectively, as shown in Fig. 7. The standard error is computed by dividing the standard deviation by the square root of the number of simulations and then representing it graphically as a shaded area. For comparison, thermal conductivity from experimental measurements, empirical potential-based MD simulation, DFT calculations, and MLIP-based MD simulations are given in Table III. Various experimental and theoretical studies have reported the κ values within the range of 84–105.5^{56,79,83–86,88} and 63–101.74 W m⁻¹ K^−157,87,88 for MoS₂ and WS₂, respectively. Remarkably, the thermal conductivity values calculated using MLIP-based MD simulations fall within this range, demonstrating favorable agreement with the DFT-calculated values. This observation underscores the capability of MLIP-based MD simulation to accurately calculate the thermal conductivity of MoS₂ and WS₂, making it a valuable tool for thermal management and precise determination of material mechanical properties, equivalent to those achieved with DFT calculations.

The interfacial thermal conductance (ITC) is the inverse of ITR. The sharp interface heterostructure was constructed to calculate ITC with a length of about 20 nm as shown in Fig. 8. Atoms in the black zone are fixed, while red and blue zones indicate the heat bath and heat sink regions, respectively. Periodic boundary condition was applied in the x (width) direction. Fixed boundary condition was applied for the direction of z. Initially, the NPT ensemble was used to equilibrate the system at 3 K and zero pressure for 30 ps with a time step of 0.1 fs. Then, the NVT ensemble was performed at 300 K and zero pressure for 800 ps with a time step of 1 fs. Next, the system was placed in the NVE ensemble for another 100 ps. After that, the hot (T_hot = 320 K) and cold (T_cold = 280 K) Nose–Hoover reservoirs were performed on the left end of the MoS₂ region and the right end of the WS₂ region for 1 ns.

The simulations were performed 8 ns under the NVE ensemble to allow the system to reach non-equilibrium steady state, where the temperature gradient along the heterostructure was well established, and the heat flux going through the system became time-independent. The temperature distribution along the heterostructure is shown in Fig. 9. The temperature in the middle region was inferred by linearly fitting the middle region of the heterostructure, and their difference was taken as ΔT. The ΔT for the interface is 13.03 K, and ITC for the MoS₂/WS₂ heterostructure is 464.89 MW m⁻² K⁻¹ calculated by Eq. (4).

In this work, we employed a deep neural network approach based on first-principles calculations to develop a universal interatomic potential tailored for modeling both thermodynamics and mechanics in MoS₂/WS₂ heterostructures. When confronted with empirical potentials in molecular simulations, the pronounced accuracy and universality of MLIPs make them significantly advantageous. This accuracy is substantiated by validating the MLIP through comparisons of RMSEs and RDFs in the MoS₂/WS₂ heterostructures against first-principles results. Using MLIP-based MD simulations, Young's modulus and lattice thermal conductivity of MoS₂ and WS₂ are predicted, aligning closely with existing literature. MLIP can be used to calculate the interfacial thermal conductance of heterostructures, while empirical potential and DFT calculations are difficult to effectively calculate it due to accuracy or computational costs. Notably, our findings indicate that the MLIP, boasting with thousands of parameters is better to empirical potentials in modeling Young's modulus and lattice thermal conductivity of MoS₂ and WS₂. Furthermore, our investigation into stress application in different directions at the interface revealed significant differences, with perpendicular stress application resulting in a notably smaller Young's modulus. This phenomenon is attributed to stress concentration at the interface. Additionally, a relative high ITC of the MoS₂/WS₂ heterostructures is found with MLIP-based NEMD simulations.

Machine learning-based atomic simulations have become a crucial tool for investigating the properties of micro- and nanostructures. Our work presents a new interatomic potential for exploring 2D TMD devices. The developed methodology also provides a valuable guide for developing MLIPs for other materials and heterostructures. Nevertheless, it is worth noting that the increased volume of parameters in MLIPs comes with higher computational costs compared to empirical potentials. In the future, enhancing the computational efficiency of MLIPs without compromising their accuracy will be a meaningful endeavor.

The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (Grant Nos. 52150610495, 12374027, and 52206080), the Shanghai Committee of Science and Technology (Grant No. 21TS1401500), and the Shanghai Pujiang Program (Grant No. 22PJ1400300).

The authors have no conflicts to disclose.

Xiangjun Liu: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal). Baolong Wang: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Kun Jia: Data curation (equal); Resources (equal); Visualization (equal). Quanjie Wang: Methodology (equal); Software (equal). Di Wang: Data curation (equal); Methodology (equal). Yucheng Xiong: Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – review & editing (equal).

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