Molecular dynamics simulations of heat transport using machine-learned potentials: A mini-review and tutorial on GPUMD with neuroevolution potentials

Heat transport properties are crucial for numerous applications.^1,2 At the atomistic level, there are primarily three computational methods for heat transport:³ molecular dynamics (MD) simulations, methods related to Boltzmann transport equation (BTE)—including, more generally, quasi-harmonic Green–Kubo (QHGK) method⁴ and Wigner transport equation (WTE) approach^5,6—combined with anharmonic lattice dynamics (ALD) (BTE-ALD for short), and atomistic Green function (AGF). Each method has its advantages and disadvantages.³ This mini-review and tutorial focuses on the MD methods. For the BTE-ALD and AGF approaches, we refer interested readers to previous tutorials.^3,7,8 Our emphasis is on thermal conductivity, including finite systems, instead of thermal boundary conductance/resistance. For the latter, we suggest referring to a previous tutorial⁹ and a review article.¹⁰ 

Notable advantages distinguish MD from the other two methods. First, MD can capture phonon–phonon scatterings at any order, while the other two methods are perturbative in nature and often consider only three-phonon scatterings (for BTE-ALD) or even completely ignore anharmonicity (for AGF). Second, MD can naturally capture scatterings of phonons by other sources such as defects and mass disorder, extending its applicability to fluid systems that are beyond the reach of the other two methods. Third, the computational cost of MD with classical potentials is usually linear with respect to the number of atoms, while it typically exhibits high-order polynomial scaling in the other two methods. Based on these considerations, MD proves particularly suitable for studying heat transport in strongly anharmonic or highly disordered systems.

Despite these advantages, MD simulations have grappled with challenges, particularly in terms of accuracy, over a considerable period of time. The predictive power of MD simulations is highly dependent on the accuracy of the classical potentials, which are mathematical models representing the potential energy surface of systems in terms of geometric information. The interatomic forces can be accurately computed using ab initio methods such as quantum-mechanical density-functional theory (DFT), leading to the ab initio molecular dynamics (AIMD) method, which has been applied to heat transport studies.^11–14 A challenge in the AIMD approach is the high computational intensity, which imposes limitations on the size and timescales that can be effectively simulated.

Recently, a type of classical potentials based on machine learning (ML) techniques, called machine-learned potentials (MLPs), has emerged as an effective framework for constructing highly accurate interatomic potentials. Due to the flexible functional forms and a large number of fitting parameters in MLPs, they can usually achieve significantly higher accuracy compared to traditional empirical potentials. Notable MLP models, to name a few, include Behler–Parrinello neural-network potential (BPNNP),¹⁵ Gaussian approximation potential (GAP) and related,^16–18 spectral neighbor analysis potential (SNAP),¹⁹ moment tensor potential (MTP),²⁰ deep potential (DP),²¹ and atomic cluster expansion (ACE).²² In this context, the recently developed neuroevolution potential (NEP) approach^23–25 simultaneously demonstrates excellent accuracy and outstanding computational efficiency, offering a distinctive advantage. Furthermore, MLPs have been increasingly used in MD simulations, including heat transport simulations (see Fig. 1 for a general trend).

Parallelization stands out as another key advancement in MD simulations, involving the deployment of parallel computing to take advantage of rapid hardware upgrades and speedups, where a large number of processors or cores work simultaneously to perform calculations, to augment computational efficiency and spatiotemporal scales of simulations. gpumd,²⁶ short for Graphics Processing Units Molecular Dynamics, represents a noteworthy development in this arena. gpumd is a versatile MD package fully implemented on graphics processing units (GPUs). This advancement facilitates the simulations of larger and more complex systems by leveraging the powerful parallel processing capabilities of GPUs. For example, it has been demonstrated that gpumd can achieve a remarkable computational speed of 1.5×10⁸ atom step/s (equivalent to a cost of $6.7× 10 − 9$ s/atom/step) in MD simulations using eight 80-gigabyte A100 graphics cards, enabling simulations up to 100 million atoms for high-entropy alloys employing a general-purpose unified NEP machine-learned potential for 16 elemental metals and their alloys.²⁷ 

In this mini-review and tutorial, we dig into the fundamentals of heat transport, the relevant MD simulation methods, and the applications of MLPs in MD simulations of heat transport. We use the NEP model^23–25 as implemented in the gpumd package²⁶ to illustrate the various technical details involved. By completing this tutorial, the readers will gain both fundamental knowledge and practical skills to construct MLPs and apply them in highly efficient and predictive MD simulations of heat transport.

When the coordinate axes align with the principal axes of the conductivity tensor, thermal transport decouples in different directions, yielding a diagonal thermal conductivity tensor with three nonzero elements: $κ x x$⁠, $κ y y$⁠, and $κ z z$⁠. These are commonly denoted as $κ x$⁠, $κ y$⁠, and $κ z$ for simplicity. For isotropic 3D systems, we usually define a conductivity scalar $κ$ in terms of the trace of the tensor: $κ=( κ x+ κ y+ κ z)/3$⁠. For isotropic 2D systems, we usually define a conductivity scalar for the planar components: $κ=( κ x+ κ y)/2$⁠. For quasi-1D systems, it is only meaningful to define the conductivity in a single direction. For simplicity, from here on we work with the conductivity scalar $κ$ unless it is necessary to consider the conductivity tensor.

At the nanoscale, the conventional concept of conductivity may lose its validity.²⁹ For example, thermal transport in materials with high thermal conductivity, such as diamond at the nanoscale, is almost ballistic, meaning the conductance changes little with increasing system length $L$⁠. In this case, if we assume that Eq. (2) still holds, then the thermal conductivity $κ$ cannot be regarded as a constant but as a function of the system length, $κ=κ(L)$⁠. This deviates from the conventional (macroscopic) concept of thermal conductivity.

We use two toy models to illustrate the above-discussed concepts. In the first model, we assume that there is only one phonon MFP of 1  $μ$m and a diffusive thermal conductivity of $κ=1000$ W m $− 1$ K $− 1$⁠. Then, ballistic conductance is $κ/λ=1$ GW m $− 2$ K $− 1$⁠. Then, the apparent thermal conductivity $κ(L)$ is given by Eq. (9), as shown in Figs. 2(a) and 2(b). In this case, $1/κ(L)$ varies linearly with $1/L$⁠. In the second model, we assume that there are two phonon modes, one with a MFP of $0.1$  $μ$m, and the other $1$  $μ$m, both having a diffusive conductivity of 500 W m $− 1$ K $− 1$⁠. Then the ballistic conductances for these two modes are 5 and 0.5 GW m $− 2$ K $− 1$⁠, respectively. The higher ballistic conductance in the second toy model can be visualized in Fig. 2(c). Although the apparent thermal conductivity for each mode follows Eq. (9), when combined, $1/κ(L)$ does not exhibit linearity with $1/L$⁠. This is an important feature for realistic materials with a general MFP spectrum $κ(ω)$⁠.

Note that the above formulation of heat current has been derived specifically for local MLPs with atom-centered descriptors. For semilocal message-passing-based MLPs, the formulation of heat current has been shown by Langer et al.^41,42 to be more complicated.

In the following, we review the heat transport MD methods implemented in the gpumd package, including equilibrium molecular dynamics (EMD), nonequilibrium molecular dynamics (NEMD), homogeneous nonequilibrium molecular dynamics (HNEMD), and spectral decomposition. While the approach-to-equilibrium method^43–45 can in principle be realized in gpumd, our discussion will primarily focus on the other three methods that have been widely employed with gpumd.

The NEMD method is a nonequilibrium and inhomogeneous method that involves implementing a pair of heat source and sink using a thermostatting method or equivalent. There are two common relative positions of the source and sink in the NEMD method, corresponding to two typical simulation setups. In one setup, the source and sink are separated by half of the simulation cell length $L$⁠, and periodic boundary conditions are applied along the transport direction. Heat flows from the source to the sink in two opposite directions in this periodic boundary setup. In the other setup, the source and sink separated by $L$ are located at the two ends of the system. Fixed boundary conditions are applied along the transport direction to prevent sublimation of the atoms in the heat source and sink. Heat flows from the source to the sink in one direction in this fixed boundary setup. It has been established⁵⁰ that the effective length in the periodic boundary setup is only $L/2$⁠. This factor must be taken into account when comparing results from the two setups.

A common practice in using the NEMD method is to extrapolate to the limit of infinite length based on the results for a few finite lengths. It is important to note that linear extrapolation is usually insufficient, as suggested even by the toy-model results shown in Fig. 2(d).

There are other spectral/modal analysis method implemented in gpumd, such as the Green–Kubo modal analysis method⁵⁹ and the homogeneous non-equilibrium modal analysis method,⁵⁸ but we will not demonstrate their usage in this tutorial.

Several MLPs have been used for heat transport with MD simulations, including BPNNP,¹⁵ GAP,¹⁶ SNAP,¹⁹ MTP,²⁰ DP,²¹ MLFF,⁶⁰ NEP,²³ SchNet,⁶¹ and So3krates.⁶²  Table I lists the relevant MLP packages implementing these MLPs.

The pioneering BPNNP model, developed by Behler and Parrinello,¹⁵ has been implemented in various packages, including runner,^15, aenet,³⁶ and kliff.⁶³ The DP, MLFF, NEP, GAP, MTP, SNAP, SchNet, and So3krates models are implemented in deepmd-kit, vasp, gpumd, quip, mlip, fitsnap, schnetpack, and mlff respectively.

Most MLP packages are interfaced to lammps³² to perform MD simulations, while NEP is native to gpumd²⁶ but can also be interfaced to lammps. The MLFF method implemented in vasp is an on-the-fly MLP that integrates seamlessly into AIMD simulations.

Table II compiles the publications up to today that have used MD simulations driven by MLPs for thermal transport studies. Note that our focus is on studies using MD simulations, excluding those solely based on the BTE-ALD approach. The number of publications up to March 10, 2024 for each MLP is shown in Fig. 1.

The application of MLPs-based MD simulations to thermal transport was pioneered by Sosso et al. in 2012 when they studied the thermal transport in the phase-changing amorphous GeTe system.⁶⁴ However, thermal transport simulations are very computationally intensive, and the rapid increase of the number of applications has only been started after the development of the GPU-based DP²¹ and NEP²³ models. In this regard, the NEP model is particularly advantageous due to its superior computational speed as compared to others.^23–25 With comparable computational resources, it has been shown to be as fast as or even faster than some empirical force fields.^154,156

There are numerous successful applications of MLPs in thermal transport. In Fig. 4, we present results from selected publications. The materials studied in these works have reliable experimental results, serving as good candidates for validating the applicability of MLPs. On one hand, MLPs demonstrate good agreement with experimental results for highly disordered materials such as liquid water,¹⁵⁴ amorphous SiO $2$⁠,¹⁴² and amorphous silicon.¹⁵² In addition to the reliability of MLPs, a crucial component for accurately describing the temperature dependence of the thermal conductivity in liquids and amorphous materials is a quantum correction method based on the spectral thermal conductivity, as defined in Eq. (39), and the quantum-statistical-correction factor, as given in Eq. (40). On the other hand, MLPs tend to systematically underestimate the thermal conductivity of crystalline solids, including silicon (using a GAP model),⁷⁴ CoSb $3$ (using a MTP model), and graphite (in-plane transport, using a NEP model).¹⁴¹ This underestimation has been attributed to the small but finite random force errors, and a correction has been devised.¹⁸¹ We will discuss this in more detail with an example in Sec. IV.

In this section, we use crystalline silicon as an example to demonstrate the workflow of constructing and using NEP models for thermal transport simulations. The NEP approach has been implemented in the open-source gpumd package.^25,26 After compiling, there will be an executable named nep that can be used to train accurate NEP models against reference data, and an executable named gpumd that can be used to perform efficient MD simulations. The gpumd package is self-contained, free from dependencies on third-party packages, particularly those related to ML. This makes the installation of gpumd straightforward and effortless. In addition, there are some handy (but not mandatory) Python packages available to facilitate the pre-processing and post-processing gpumd inputs and outputs, including calorine,^182, gpyumd,^183, gpumd-wizard,¹⁸⁴ and pynep.¹⁸⁵ Since its inception with the very first version in 2013,¹⁸⁶  gpumd has been developed with special expertise in heat transport applications.

The NEP model is based on artificial neural network (ANN) and is trained using a separable natural evolution strategy (SNES),¹⁸⁷ hence the name.

The free parameters are optimized using the SNES by minimizing a loss function that is a weighted sum of the root-mean-square errors (RMSEs) of energy, force, and virial stress, over $N gen$ generations with a population size of $N pop$⁠. The weights for the energy, force, and virial are denoted $λ e$⁠, $λ f$⁠, and $λ v$⁠, respectively. Additionally, there are proper norm-1 (⁠ $ℓ 1$⁠) and norm-2 (⁠ $ℓ 2$⁠) regularization terms. For explicit details on the training algorithm, refer to Ref. 23.

Although NEP with proper hyperparameters can account for almost all types of interactions, it can be useful to combine it with some well developed potentials, such as the Ziegler–Biersack–Littmark (ZBL)¹⁸⁸ potential for describing the extremely large screened nuclear repulsion at short interatomic distances and the D3 dispersion correction¹⁸⁹ for describing relatively long-range but weak interactions. Both potentials have been recently added to the gpumd package.^168,190 It has been demonstrated that dispersion interactions can reduce the thermal conductivity of typical metal-organic frameworks by about 10%.¹⁶⁸ With the addition of ZBL and D3, NEP can then focus on describing the medium-range interactions.

There are educational articles focusing on various best practices in constructing MLPs.^191,192 Here, we use crystalline silicon as a specific example to illustrate the particular techniques in the context of NEP.

A training dataset is a collection of structures, each characterized by a set of attributes:

1. a cell matrix defining a periodic domain

2. the species of the atoms in the cell

3. the positions of the atoms

4. the total energy of the cell of atoms

5. the force vector acting on each of the atoms

6. (optionally) the total virial tensor (with six independent components) of the cell

While there are already several publicly available training datasets for silicon, we opt to create one from scratch for pedagogical purposes. The construction of training dataset typically involves an iterative process, employing a scheme similar to active learning. The iterative process begins with an initial dataset. To investigate heat transport in crystalline silicon, the initial training dataset should encompass structures relevant to the target temperatures and pressures. The most reliable way of generating structures under these conditions is through performing AIMD simulations, where interatomic forces are calculated based on quantum mechanical methods, such as the DFT approach. However, AIMD is computationally expensive (which is the primary motivation for developing a MLP) and it is often impractical to perform AIMD simulations for a dense grid of thermodynamic conditions. Fortunately, there is usually no such need for the purpose of generating the reference structures. Actually, manual perturbation of the atomic positions and/or the cell matrices proves to be an effective way of generating useful reference structures.

Based on the considerations above, we generate the initial training dataset through the following methods. First, we generate 50 structures by applying random strains (ranging from $−3%$ to $+3%$ for each degree of freedom) to the unit cell of cubic silicon (containing 8 atoms) while simultaneously perturbing the atomic positions randomly (by 0.1 Å). Second, we perform a 10-ps AIMD simulation at 1000 K (fixed cell) using a $2×2×2$ supercell of silicon crystal containing 64 atoms, and sample the structures every 0.1 ps, obtaining another 100 structures. In total, we obtain 150 structures and 6800 atoms initially.

After obtaining the structures, we perform single-point DFT calculations to obtain the reference energy, force, and virial data. These data are saved to a file named train.xyz, using the extended XYZ format. The single-point DFT calculations are performed using the vasp package,¹⁹³ using the Perdew–Burke–Ernzerhof functional with the generalized gradient approximation,¹⁹⁴ a cutoff energy of 600 eV, an energy convergence threshold of $10 − 6$ eV, and a $k$-point mesh of $4×4×4$ for 64-atom supercells and $12×12×12$ for 8-atom unit cells.

With the training data, we proceed to train our first NEP model, denoted as NEP-iteration-1. For this task, we need to prepare an input file named nep.in for the nep executable in the gpumd package. This nep.in input file contains the various hyperparameters for the NEP model under training. Most hyperparameters have well-suited default values, and for users initiating this process, it is recommended to use these defaults whenever applicable. The default values for key hyperparameters are as follows:

1. $n max R:4$

2. $n max A:4$

3. Chebyshev polynomial basis size for radial descriptor components $N bas R:12$

4. Chebyshev polynomial basis size for angular descriptor components $N bas A:12$

5. $l max 3 body:4$

6. $l max 4 body:2$

7. $l max 5 body:0$ (not used by default)

8. $N neu:30$

9. Energy and force weights $λ e$ and $λ f:1$

10. Virial weight $λ v:0.1$

11. Batch size: 1000 (a large or even full batch is preferred for training with SNES)

12. Population size in SNES: 50

13. Number of training generations (steps): 100 000

14. ANN architecture: 30-30-1 (input layer 30, hidden layer 30, scalar output; relatively small but sufficient for most cases, expect for very complicated training data.)

Following this strategy, we use a very simple nep.in input file for our case, which is as follows:

The training results for NEP-iteration-1 are shown in Fig. 5(a). The RMSEs of force, energy, and virial all converge well within the default 100 000 training steps. The parity plots for force, energy, and virial in Figs. 5(b)–5(d) show good correlations between the NEP predictions and the DFT reference data. The RMSEs for energy, force, and virial are 1.0 meV/atom, 54.6 meV/Å, and 21.8 meV/atom, respectively.

Reliable assessment of the accuracy of a MLP typically involves an independent test dataset rather than the training dataset. To this end, we perform 10-ps MD simulations using NEP-iteration-1 in the $NPT$ ensemble. The target pressure is set to zero, and the target temperatures range from 100 to 1000 K with intervals of 100 K. We sample 100 structures, totalling 6400 atoms.

We perform single-point DFT calculations for these structures and then use NEP-iteration-1 to generate predictions. This is achieved by adding the prediciton keyword to the nep.in file:

The high accuracy of the latest test dataset sampled from 10-ns MD simulations driven by NEP-iteration-2 suggests that NEP-iteration-2 is a reliable model for MD simulation of crystalline silicon from 100 to 1000 K. Therefore, we conclude the iteration and use NEP-iteration-2 for the thermal transport applications. In the following, we will refer to NEP-iteration-2 simply as NEP. This NEP model, running on a consumer-grade NVIDIA RTX 4090 GPU card with 24 GB of memory, achieves a remarkable computational speed of about $2.4× 10 7$ atom-step/second, equivalent to a computational cost of about $4.2× 10 − 8$ s/atom/step in MD simulations.

Using a trained MLP to generate MD trajectory is a common practice in nearly all the active-learning schemes documented in the literature. The major difference between different active-learning schemes is about the criteria for selecting structures to be added to the training dataset. While there might be a risk of sampling nonphysical structures using a trained MLP model, as demonstrated in this tutorial, one can mitigate the risk by conducting a few iterations and employing shorter MD runs in the initial stages, progressively increasing the MD simulation time with each iteration. As a result, the MLP becomes increasingly reliable throughout the iteration process, enabling the generation of longer and more accurate trajectories over time. In our example using the silicon crystal, a relatively simple system, we only performed two iterations to achieve accurate predictions for 10-ns MD runs. However, for more complex systems, one might need to perform more iterations, increasing the MD steps more gradually than what we have demonstrated for the silicon crystal example.

Before applying a MLP to thermal transport applications, it is usually a good practice to examine the phonon dispersion relations. The phonon dispersion relations for NEP and Tersoff¹⁹⁵ potentials are calculated using gpumd, employing the finite-displacement method with a displacement of 0.01 Å. For DFT, we use density functional perturbation theory as implemented in vasp in combination with phonopy,¹⁹⁶ using a $4×4×4$ supercell, a cutoff energy of 600 eV, an energy convergence threshold of $10 − 8$ eV, and a $5×5×5$ $k$-point mesh.

In Fig. 6, we compare the phonon dispersion relations calculated from DFT, Tersoff potential, and NEP. While there are small differences between NEP and DFT results, the agreement between NEP and DFT is significantly better than that between Tersoff and DFT. The agreement between NEP and DFT can, in principle, be further improved, for example, by increasing the size of the ANN model and/or the cutoff radii. However, this comes with a trade-off, as it may reduce computational efficiency. In practice, achieving a balance between accuracy and speed is essential.

After validating the phonon dispersion relations, we proceed to thermal conductivity calculations using the various MD methods, as reviewed in Sec. III. All calculations are performed using the gpumd executable in gpumd.

We start with the EMD method, using a sufficiently large $12×12×12$ cubic supercell with 13 824 atoms. The run.in file for the gpumd executable is configured as follows:

We perform 50 independent runs using the inputs above, each with a different set of initial velocities. The $κ(t)$ [cf. Eq. (26)] results from individual runs (thin solid lines) and their average (thick solid line) and error bounds (thick dashed lines) are shown in Fig. 7(a). Taking $t=1$ ns as the upper limit of the correlation time, up to which $κ(t)$ converges well, we have $κ≈102±6$ W m $− 1$K $− 1$ from the EMD method. In this work, all statistical errors are calculated as the standard error of the mean.

We then move to the HNEMD method. Since the HNEMD method has the same finite-size effects as in the EMD method, we use the same simulation cell as in the EMD method. The run.in file for the gpumd executable reads as follows:

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.⁵⁴ 

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)].

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.

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.

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.

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.

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).

The authors have no conflicts to disclose.

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).

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.

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.