We perform 4 independent runs using the specified inputs, each with a different set of initial velocities. The $κ ( t )$ [cf. Eq. (34) ] results from individual runs (thin solid lines) and their average (thick solid line) and error bounds (thick dashed lines) are shown in Fig. 7(b) . The estimated thermal conductivity is $κ ≈ 108 ± 4$ Wm $− 1$ K $− 1$ ⁠ , consistent with the EMD value within statistical error bounds. It is noteworthy that the total production time for the HNEMD simulations ( ⁠ $4 × 10$ ns) is considerably smaller than that for the EMD simulations ( ⁠ $50 × 10$ ns), while the former still gives a smaller statistical error. This suggests a higher computational efficiency of the HNEMD over the EMD method, as previously emphasized. ⁵⁴

4. Thermal conductivity from NEMD

The third MD method we demonstrate is the NEMD method, using the fixed boundary setup discussed in Sec. II C 2 . We explore lengths $L = 2.7$ ⁠ , 5.5, 11.0, 21.9, 43.8, 87.6, 175.3, 350.5 nm, maintaining a consistent $5 × 5$ cell in the transverse direction. The heat source and sink regions span 4.4 nm, which is long enough to ensure fully thermalized phonons within these regions. The run.in input file for our NEMD simulations reads as follows:

Unlike the EMD and HNEMD simulations, the NEMD simulations involve an extra operation: certain atoms are frozen. We assign these atoms to the “group” 0 and use the fix 0 command to freeze them. In the production stage, two Langevin thermostats with different temperatures are applied separately to groups 1 and 7, corresponding to the heat source and the heat sink, respectively. The temperature difference between them is set to 20 K. The heat flux can be obtained from the data produced by the compute keyword, allowing us to calculate the apparent thermal conductivity $κ ( L )$ according to Eq. (29) . The production stage has a duration of 2 ns, with a well-established steady state achieved within the first 1 ns. Therefore, we use the second half of the production time to calculated the aforementioned stead-state properties. For each system length, we perform 2 independent runs, each with a different set of initial velocities. To get the spectral conductance $G ( ω )$ in the ballistic limit, as used in Eq. (45) , we use the data produced by the compute_shc keyword in NEMD simulations with a short system length of $L = 1.6$ nm.

As expected, the $κ ( L )$ values from NEMD simulations match well with the $κ ( L )$ curve from the HNEMD-based formalism [ Fig. 7(d) ]. However, reaching the diffusive limit directly through NEMD simulations is computationally demanding. Considering the presence of different phonon MFPs [ Fig. 7(c) ] in the system, linear extrapolation to the diffusive limit based on a limited number of $κ ( L )$ values from NEMD simulations is often inadequate. This limitation arises because the relation between $1 / κ ( L )$ and $1 / L$ becomes nonlinear in the large- $L$ limit (see Fig. 8 ). This nonlinearity is a general feature in realistic materials, as also demonstrated in our toy model [ Fig. 2(d) ].

FIG. 8.

View large Download slide

The nonlinearity in the relation between $κ ( L = ∞ ) / κ ( L )$ and $1 / L$ in the large- $L$ limit, observed in the second toy model [as discussed in Fig. 2(d) ] and the silicon example.

FIG. 8.

View large Download slide

The nonlinearity in the relation between $κ ( L = ∞ ) / κ ( L )$ and $1 / L$ in the large- $L$ limit, observed in the second toy model [as discussed in Fig. 2(d) ] and the silicon example.

Close modal

As of now, we have demonstrated the full consistency among the three MD-based methods. Notably, the HNEMD method stands out as the most computationally efficient. This explains why most works based on gpumd utilize the HNEMD method, with the other two methods typically being employed primarily for sanity-checking the results.

5. Comparison with experiments

After obtaining consistent results from three MD methods, we are ready to compare the results with experimental data. The thermal conductivity of crystalline silicon is measured to be about 150 W m $− 1$ K $− 1$ ⁠ , but our HNEMD simulations predict a value of $108 ± 4$ W m $− 1$ K $− 1$ ⁠ , which is only $72 %$ of the experimental value. As a comparison, the thermal conductivity of crystalline silicon has been calculated ¹⁹⁷ to be about $250 ± 10$ W m $− 1$ K $− 1$ using a Tersoff potential, ¹⁹⁵ which is $167 %$ of the experimental value. Specifically, the Tesoff potential appears to underestimate the phonon anharmonicity, while the NEP model tends to overestimate it.

Based on the correction method, we perform HNEMD simulations using the Langevin thermostat with the following set of $τ T$ values: 30, 50, 100, 200, and 500 ps. From these, the $σ L$ values are calculated to be 17.3, 27.4 38.7 54.8, and 70.7 meV/Å. At 300 K, the force RMSE for our NEP model is tested to be $σ mlp = 21.2$ meV/Å. Therefore, the resulting $σ tot$ values are 27.4, 34.6, 44.2, 58.7, and 73.8 meV/Å. To ensure consistency with experimental conditions, we also account for the presence of a few Si isotopes (92.2% $28$ Si, 4.7% $29$ Si, and 3.1% $30$ Si) in the calculations. The calculated $κ ( σ tot )$ with the various $σ tot$ are shown in Fig. 9(a) . By fitting these data, we obtain a corrected thermal conductivity of $κ ( σ tot = 0 ) = 145$ W m $− 1$ K $− 1$ ⁠ , in excellent agreement with the experimental value.

FIG. 9.

View large Download slide

(a) Thermal conductivity of crystalline silicon at 300 K from HNEMD simulations using the herein developed NEP models as a function of the total force error $σ tot$ ⁠ . NHC and LAN represent the Nosé–Hoover chain and Langevin thermostatting methods, respectively. The data are fitted to obtain the corrected thermal conductivity of $κ ( σ tot = 0 ) = 145$ W m $− 1$ K $− 1$ ⁠ . (b) Comparison of simulation results before and after the correction with experimental values ^178,198,199 and previously (uncorrected and corrected) NEP-MD simulations. ¹⁸¹

FIG. 9.

View large Download slide

(a) Thermal conductivity of crystalline silicon at 300 K from HNEMD simulations using the herein developed NEP models as a function of the total force error $σ tot$ ⁠ . NHC and LAN represent the Nosé–Hoover chain and Langevin thermostatting methods, respectively. The data are fitted to obtain the corrected thermal conductivity of $κ ( σ tot = 0 ) = 145$ W m $− 1$ K $− 1$ ⁠ . (b) Comparison of simulation results before and after the correction with experimental values ^178,198,199 and previously (uncorrected and corrected) NEP-MD simulations. ¹⁸¹

Close modal

The extrapolation scheme described by Eq. (48) not only applies to a single NEP model with different levels of intentionally added random forces through the Langevin thermostat, but is also valid for different NEP models with varying accuracy. To demonstrate this, we construct two extra NEP models with reduced accuracy. Starting from the default hyperparameters, we construct the first extra NEP model by reducing the number of neurons in the hidden layer from 30 to 1, resulting in an increased force RMSE of 32.4 meV/Å. Based on this, we then construct the second extra NEP model by further reducing the Chebyshev polynomial basis sizes ( ⁠ $N bas R$ ⁠ , $N bas A$ ⁠ ) from $( 12 , 12 )$ to $( 4 , 4 )$ ⁠ , resulting in a further increased force RMSE of 52.9 meV/Å. The thermal conductivity results from the three NEP models with different accuracy using the Nosé–Hoover chain thermostat also closely follow the extrapolation curve [ Fig. 9(a) ], providing further support for the validity of the extrapolation scheme Eq. (48) .

Our results for 300 K before and after the correction are consistent with those reported in the previous work, ¹⁸¹ which also uses a NEP model [ Fig. 9(b) ]. In Fig. 9(b) , we also show the results for other temperatures ¹⁸¹ in comparison to the experimental data. The corrected results agree well with the experimental ones across a broad range of temperatures. The slightly higher values from corrected NEP model predictions are likely due to the fact that isotope disorder was not considered in the previous calculations. ¹⁸¹

While we have only demonstrated the application of the extrapolation (correction) method to HNEMD simulations, it is worth noting that this method is also potentially applicable to EMD simulations. We speculate that the force errors in MLPs may also play a role in ALD-based approaches for thermal transport.

V. SUMMARY AND CONCLUSIONS

In summary, we have provided a comprehensive pedagogical introduction to MD simulations of thermal transport utilizing the NEP MLP as implemented in the gpumd package.

We began by reviewing fundamental concepts related to thermal transport in both ballistic and diffusive regimes, elucidating the explicit expression of the heat flux in the context of MLPs, and exploring various MD-based methods for thermal transport studies, including EMD, NEMD, HNEMD, and spectral decomposition.

Following this, we conducted an up-to-date review of the literature on the application of MLPs in thermal transport problems through MD simulations.

A detailed review of the NEP approach followed, with a step-by-step demonstration of the process of developing an accurate and efficient NEP model for crystalline silicon applicable across a range of temperatures. Utilizing the developed NEP model, we explained the technical details of all MD-based methods for thermal transport discussed in this work. Finally, we compared the simulation results with experimental data, addressing the common trend of thermal conductivity underestimation by MLPs and demonstrating an effective correction method.

By completing this tutorial, readers will be equipped to construct MLPs and seamlessly integrate them into highly efficient and predictive MD simulations of heat transport.

ACKNOWLEDGMENTS

H.D. is supported by the Science Foundation from the Education Department of Liaoning Province (No. JYTMS20231613) and the Doctoral start-up Fund of Bohai University (No. 0523bs008). P.Y. is supported by the Israel Academy of Sciences and Humanities & Council for Higher Education Excellence Fellowship Program for International Postdoctoral Researchers. K.X. and T.L. acknowledge support from the National Key R&D Project from Ministry of Science and Technology of China (No. 2022YFA1203100), the Research Grants Council of Hong Kong (No. AoE/P-701/20), and RGC GRF (No. 14220022). Z.Z. acknowledges the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 101034413. S.X. acknowledges financial support from the National Natural Science Foundation of China (NNSFC) (Grant No. 12174276).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Haikuan Dong: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Yongbo Shi: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Penghua Ying: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Ke Xu: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Ting Liang: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Yanzhou Wang: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Zezhu Zeng: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal). Xin Wu: Conceptualization (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Wenjiang Zhou: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal). Shiyun Xiong: Conceptualization (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Shunda Chen: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Zheyong Fan: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

All the training and test datasets and the trained NEP models for crystalline silicon are freely available at https://gitlab.com/brucefan1983/nep-data . The training datasets, trained NEP, DP, and MTP models for graphene and MD input files for reproducing Fig. 3 are freely available at https://github.com/hityingph/supporting-info/tree/main/Dong_GPUMD_Tutorial_2024 .

APPENDIX: DETAILS ON HEAT CURRENT VALIDATION

To validate the implementation of heat current for NEP, DP, and MTP (see Fig. 3 ), we use a common reference dataset to train a model for each of the three MLPs. We take all the monolayer graphene structures from Ref. 67 and use 3288 structures (99 493 atoms) as our training dataset and 822 structures (25 035 atoms) as our test dataset, respectively.

For NEP, we set $r c R = 6$ Å, $r c A = 4$ Å, $N neu = 50$ ⁠ , $N gen = 5 × 10 5$ ⁠ , while keeping other hyperparameters as the defaults. For DP, the deepmd-kit package (version 2.1.4) ²¹ is used, with the se_a descriptor with a cutoff radius of 6 Å. The dimensions of the embedding network are set to (25, 50, 100), and the fitting network dimensions are configured as (240, 240, 240). Initially, the weighting parameters for energy and forces are set to 0.02 and 1000, respectively, and are linearly adjusted to 1 for both during the training process. The training comprises $4 × 10 6$ steps, with a learning rate that is exponentially decreased from $10 − 3$ to $10 − 8$ ⁠ . For MTP, the mlip (version 2) package ²⁰ is used. The descriptor “level” for MTP is set to 18, with a cutoff radius of 6 Å. Table III presents the performance metrics for the three MLP models.

We then conduct NEMD simulations to validate the implementations of heat current in the three MLPs by checking the consistency between the accumulated heat in atoms within the transport region [cf. Eq. (25) ] and that obtained from the thermostats [cf. Eq. (31) ]. The NEMD simulation procedure is similar to that as described in Sec. IV C 4 for silicon. The transport is set to be along the armchair direction of a graphene sample with a width of 2.5 nm and a length of 426 nm (excluding the thermostatted regions). The data presented in Fig. 3 are sampled during the last 1.5 ns of the NEMD simulations, during which a steady state is achieved.

REFERENCES