Molecular

Molecular dynamics (MD) simulations play an important role in understanding and engineering heat transport properties of complex materials. An essential requirement for reliably predicting heat transport properties is the use of accurate and efficient interatomic potentials. Recently, machine-learned potentials (MLPs) have shown great promise in providing the required accuracy for a broad range of materials. In this mini review and tutorial, we delve into the fundamentals of heat transport, explore pertinent MD simulation methods, and survey the applications of MLPs in MD simulations of heat transport. Furthermore, we provide a step-by-step tutorial on developing MLPs for highly efficient and predictive heat transport simulations, utilizing the neuroevolution potentials (NEPs) as implemented in the GPUMD package. Our aim with this mini review and tutorial is to empower researchers with valuable insights into cutting-edge methodologies that can significantly enhance the accuracy and efficiency of MD simulations for heat transport studies.


I. INTRODUCTION
Heat transport properties are crucial for numerous applications 1,2 .
At the atomistic level, there are primarily three computational methods for heat transport 3 : molecular dynamics (MD) simulations, methods related to Boltzmann transport equation (BTE)-including, more generally, quasi-harmonic Green-Kubo (QHGK) method 4 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 a) Electronic mail: donghaikuan@163.comb) Electronic mail: phychensd@gmail.com c Electronic mail: brucenju@gmail.comand disadvantages 3 .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 9 and a review article 10 .
Notable advantages distinguish MD from the other two methods.Firstly, 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 the anharmonicity (for AGF).Secondly, 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.Thirdly, 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][12][13][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) 15 , Gaussian approximation potential (GAP) and related [16][17][18] , spectral neighbor analysis potential (SNAP) 19 , moment tensor potential (MTP) 20 , deep potential (DP) 21 , and atomic cluster expansion (ACE) 22 .In this context, the recently developed neuroevolution potential (NEP) approach [23][24][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 26 , 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 8 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 highentropy alloys employing a general-purpose unified NEP machine-learned potential for 16 elemental metals and their alloys 27 .
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][24][25] as implemented in the gpumd package 26 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.

Thermal conductivity
Fourier's law describes the empirical relationship governing heat transport, expressed as: Here Q µ is the heat flux in the µ direction, ∂T ∂xν is the temperature gradient in the ν direction, and κ µν is an element of the second-rank conductivity tensor 28 .The heat flux measures the heat transport per unit time and per unit area, typically measured in W m −2 .The thermal conductivity is commonly expressed in units of W m −1 K −1 .
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: κ xx , κ yy , and κ zz .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.

Thermal conductance
In macroscopic transport (the meaning of which will become clear soon), thermal conductance K is related to thermal conductivity by: where A is the cross-sectional area and L is the system length along the transport direction.This relation is similar to that between electrical conductance and electrical conductivity one learns in high school.Usually, thermal conductivity is considered an intrinsic property of a material, and thermal conductance depends on the geometry (A and L).However, complexities emerge when examining heat transport at the nanoscale or mesoscale.At the nanoscale, the conventional concept of conductivity may lose its validity 29 .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.
Rather than adhering strictly to Eq. ( 2), one can generalize the relation between conductance and conductivity as follows: where K 0 is the ballistic thermal conductance of the material.The term κ in Eq. (3) refers to the diffusive thermal conductivity, the conventional thermal conductivity defined in the macroscopic limit (L → ∞) where the phonon transport is diffusive.By contrast, the lengthdependent thermal conductivity κ(L) defined in Eq. ( 2) is usually called the apparent thermal conductivity or effective thermal conductivity.In the diffusive limit, the apparent thermal conductivity κ(L) defined in Eq. ( 2) approaches the diffusive conductivity κ defined in Eq.
(3), as expected.By comparing Eq. ( 2) and Eq. ( 3), we obtain the following relation between the apparent thermal conductivity κ(L) and the diffusive thermal conductivity κ: From this, we have It is more common to use thermal conductance per unit area G, which is defined as The corresponding ballistic conductance per unit area is The ratio between the diffusive conductivity and the ballistic conductance per unit area defines a phonon mean free path (MFP): In terms of the phonon MFP, we have This is known as the ballistic-to-diffusive transition formula for the length dependent thermal conductivity.Figure 2 schematically shows the ballistic-to-diffusive transition behavior.The above discussion is simplified in the sense that no channel dependence of the thermal transport has been taken into account.Different channels usually have different MFPs and diffusive conductivities.In general, both the conductivity and the MFP are frequency dependent and we can generalize Eq. ( 9) to With κ(ω, L), we can obtain the apparent thermal conductivity at any length L as 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 the 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 GW m −2 K −1 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 κ(ω).

B. Heat flux and heat current
The heat flux is defined as the time derivative of the sum of the moments of the site energies of the particles in the system 30 : The site energy E i is the sum of the kinetic energy m i v 2 i /2 and the potential energy U i .Here m i , r i , and v i are the mass, position, and velocity of particle i, respectively, and V is the controlling volume for the particles, which is usually the volume of the simulation box, but can also be specifically defined for low-dimensional systems simulated with vacuum layers.In MD simulations, it is usually more convenient to work on the heat current that is an extensive quantity: It is clear that the total heat current can be written as two terms: where the first term is the kinetic or convective part, and the second term is called the potential part, The expression above involves absolute positions and is thus not directly applicable to periodic systems.To derive an expression that can be used for periodic systems, we need to discuss potential energy and interatomic force first.
For the MLPs discussed in this tutorial, the total potential energy U of a system can be written as the sum of site potentials U i : The site potential can have different forms in different potential models.A well-defined force expression for general many-body potentials that explicitly respects Newton's third law has been derived as 31 : where Here, ∂U i /∂r ij is a shorthand notation for a vector with Cartesian components ∂U i /∂x ij , ∂U i /∂y ij , and ∂U i /∂z ij .The atomic position difference is defined as Using the force expression, the heat current can be derived to be 31 : From the definition of virial tensor and the force expression Eq. ( 18), we have Using the explicit force expression Eq. ( 19), we can also express the per-atom virial as Therefore, the heat current can be neatly written as: This expression, which involves relative atom positions only, is applicable to periodic systems and has been implemented in the gpumd package 26 for all the supported interatomic potentials, including NEP.The current implementation of the heat current in lammps 32 is generally incorrect for many-body potentials, and corrections to lammps have only been done for special force fields 33,34 .For any MLP that interfaces with lammps, one must use the full 9 components of the per-atom virial and provide a correct implementation of Eq. ( 24).NEP has an interface for lammps that meets this requirement.To the best of our knowledge, among the other publicly available MLP packages, only deepmd 21 (after the work of Tisi et al. 35 ) and aenet 36 (after the work of Shimamura et al. 37 ) have implemented the heat current correctly.The heat current is also correctly formulated 38 for a MLP based on the smooth overlap of atomic positions 39 .Contrarily, the widely used MTP method 20 (as implemented in Ref. 40), for example, exhibits an incorrect implementation of the heat current, as demonstrated in Fig. 3.According to energy conservation, the accumulated heat from the atoms [cf.Eq. ( 25)] should match that from the thermostats [cf.Eq. ( 31)], allowing for only small fluctuations.It is evident that both DP and NEP exhibit this property, whereas MTP does not.Details on the calculations are provided in Appendix A. Note that the above formulation of heat current has been derived specifically for local MLPs with atomcentered descriptors.For semilocal message-passingbased MLPs, the formulation of heat current has been shown by Langer et al. 41,42 to be more complicated.

C. Overview of MD-based methods for heat transport
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][44][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 EMD method
The EMD method is based on the Green-Kubo relation for thermal transport 46 : where C µν (t) is the heat current autocorrelation function (HCACF) The equations above define the running thermal conductivity, which is a function of the correlation time t.In MD simulations, the correlation function is defined as where t p is the production time within which the heat current data are sampled.This production run should be in an equilibrium ensemble (as indicated by the subscript "e" in the HCACF expression), usually N V E, but N V T with a global thermostat can also be used.Thermal conductivity in the diffusive limit is obtained by taking the limit of t → ∞, but in practice, this limit can be well approximated at an appropriate t.One also needs to ensure that the simulation cell is sufficiently large to eliminate finite-size effects 47-49 .

The NEMD method
The NEMD method is a nonequilibrum 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 50 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.
When the system reaches a steady state, a temperature profile with a definite temperature gradient ∇T will be established.Meanwhile, a steady heat flux Q will be generated.With these, one can obtain the apparent thermal conductivity κ(L) of a system of finite length L according to Fourier's law, in the linear response regime where the temperature gradient |∇T | across the system is sufficiently small.It has been observed that the local Langevin thermostat outperforms the global Nosé-Hoover thermostat 51,52 in generating temperature gradients 53 .It has also been demonstrated that the temperature gradient should be directly calculated from the temperature difference |∇T | = ∆T /L rather than through fitting part of the temperature profile 53 .This is to ensure that the contact resistance is also included, and the total thermal conductance is given by The steady-state heat flux can be computed either microscopically or from the energy exchange rate dE/dt in the thermostatted regions and cross-sectional area A as based on energy conservation.The two approaches must generate the same result, and they have been used to validate the implementation of heat flux in several MLPs, as shown in Fig. 3.
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).

The HNEMD method
In the HNEMD method, an external force of the form 54 is added to each atom to drive the system out of equilibrium, inducing a nonequilibrum heat current (note the subscript "ne"): The driving force parameter F e is of the dimension of inverse length.The quantity in the parentheses is proportional to the running thermal conductivity tensor and we have This provides a way of computing the thermal conductivity.In the HNEMD method, the system is in a homogeneous nonequilibrium state because there is no explicit heat source and sink.The system is periodic in the transport direction and heat flows circularly under the driving force.Because of the absence of heat source and sink, no boundary scattering occurs for the phonons and the HNEMD method is similar to the EMD method in terms of finite-size effects.

Spectral decomposition
In the framework of the NEMD and HNEMD methods, one can also calculate spectrally decomposed thermal conductivity (or conductance) using the virial-velocity correlation function 54,55 In terms of this, the thermal conductance in NEMD simulation can be decomposed as follows: The thermal conductivity in HNEMD simulation can be decomposed as follows: The virial-velocity correlation function here is essentially the force-velocity correlation function defined for a (physical or imaginary) interface 56,57 .The spectral quantities allow for a feasible quantumstatistical correction 3,58 for strongly disordered systems where phonon-phonon scatterings are not dominant.For example, the spectral thermal conductivity can be quantum-corrected by multiplying the factor where x = ℏω/k B T .There are other spectral/modal analysis method implemented in gpumd, such as the Green-Kubo modal analysis method 59 and the Homogeneous non-equilibrium modal analysis method 58 , but we will not demonstrate their usage in this tutorial.

III. REVIEW OF MD SIMULATION OF HEAT TRANSPORT USING MLPS
Several MLPs have been used for heat transport with MD simulations, including BPNNP 15 , GAP 16 , SNAP 19 , MTP 20  Most MLP packages are interfaced to lammps 32 to perform MD simulations, while NEP is native to gpumd 26 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 10th, 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 phasechanging amorphous GeTe system 64 .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 GPUbased DP 21 and NEP 23 models.In this regard, the NEP model is particularly advantageous due to its superior computational speed as compared to others [23][24][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 154 , amorphous SiO 2 142 , and amorphous silicon 152 .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) 74 , CoSb 3 (using a MTP model), and graphite (inplane transport, using a NEP model) 141 .This underestimation has been attributed to the small but finite random force errors, and a correction has been devised 181 .We will discuss this in more detail with an example in the next section.

IV. MOLECULAR DYNAMICS SIMULATION OF HEAT TRANSPORT USING NEP AND GPUMD
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 TABLE I.The MLPs and their implementation packages that have been used in MD simulations of heat transport.

Year MLP Package
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 selfcontained, 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 184 , and pynep 185 .Since its inception with the very first version in 2013 186 , gpumd has been developed with special expertise in heat transport applications.

A. The neuroevolution potential
The NEP model is based on artificial neural network (ANN) and is trained using a separable natural evolution strategy (SNES) 187 , hence the name.

The NN model
The ML model in NEP is a fully-connected feedforward ANN with a single hidden layer, which is also called a multilayer perceptron.The total energy is the sum of the site energies U = i U i , and the site energy U i is the output of the neural network (NN), expressed as: ω (1) µ tanh Here, N des is the number of descriptor components, N neu is the number of neurons in the hidden layer, q i ν is the ν-th descriptor component of atom i, ω (0) µν is the connection weight matrix from the input layer to the hidden layer, ω µ is the connection weight vector from the hidden layer to the output layer, b (0) µ is the bias vector in the hidden layer, and b (1) is the bias in the output layer.ω µ , and b (1) are trainable parameters.The function tanh(x) is the nonlinear activation function in the hidden layer.According to Eq. ( 41), the NEP model is a simple analytical function of a descriptor vector.A C++ function for evaluating the energy and its derivative with respect to the descriptor components can be found in Ref. 25.

The descriptor
The descriptor q ν i encompasses the local environment of atom i.In NEP, the descriptor is an abstract vector whose components group into radial and angular parts.The radial descriptor components q i n (0 ≤ n ≤ n R max ) are defined as where r ij is the distance between atoms i and j and g n (r ij ) are a set of radial functions, each of which is formed by a linear combination of Chebyshev polynomials.The angular components include n-body (n = 3, 4, 5) correlations.For the 3-body part, the descriptor components are defined as (0 Here, Y lm are the spherical harmonics and rij is the unit vector of r ij .Note that the radial functions g n (r ij ) for the radial and angular descriptor components can have different cutoff radii, which are denoted as r R c and r A c , respectively.For 4-body and 5-body descriptor components (with similar hyperparameters l 4body  Diamond allotropes 2023 Ying 156,157 MOF crystals, Phosphorene Zhang 158 Amorphous HfO2 2024 Cao 159 Phosphorous carbide Cheng 160 Perovskites 2024 Fan 161 HKUST-1 crystal Fan 162 Graphene antidot lattice 2024 Li 163 Strained monolayer graphene Li 164 Amorphous silicon 2024 Li 165 2D COF-5 Pegolo 166 Glassy LixSi1−x 2024 Wang 167 Ga2O3 Ying 168 MOF crystals 2024 Yue 169 Si-C interfaces Zeraati 170 La2Zr2O7 and many others 2024 Zhang 171 GeTe So3krates Langer 42 SnSe

The training algorithm
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.

Combining with other potentials
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) 188 potential for describing the extremely large screened nuclear repulsion at short interatomic distances and the D3 dispersion correction 189 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% 168 .With the addition of ZBL and D3, NEP can then focus on describing the medium-range interactions.

B. Model training and testing
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.

Prepare the initial training data
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 6 independent components) of the cell The structures can be prepared by any method, while the energy, force, and virial are usually calculated via quantum mechanical methods, such as the DFT method.For a dataset comprising N str structures with a total number of N atoms, there are N str energy data, 6N str virial data, and 3N force data.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.Firstly, 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 Å).Secondly, 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 singlepoint 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 193 , using the Perdew-Burke-Ernzerhof functional with the generalized gradient approximation 194 , 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.

Train the first NEP model
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.infor 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: Following this strategy, we use a very simple nep.in input file for our case, which is as follows: In the first line, we specify the number of species (atom types) and the chemical symbol(s).In our example, there is only one species with the chemical symbol Si.In the second line, we specify the cutoff radii r R c and r A c for the g n (r ij ) functions in the radial and angular descriptor components, respectively.In our example, both cutoff radii are set to 5 Å, which includes the third nearest neighbors.The choice of cutoff radii is crucial for the performance of the trained NEP model and usually requires a systematic exploration to find an optimal set of values.It is important to note that the average number of neighbors, and hence the computational cost, scales cubically with respect to the cutoff radii.Therefore, blindly using large cutoff radii is not advisable.Although r R c = r A c in our current example, it is generally beneficial to use a larger r R c and a smaller r A c , because the radial descriptor components are computationally much cheaper than the angular descriptor components.Using a larger r R c does not lead to a significant increase in the computational cost, but can help capture longer-range interactions (such as screened Coulomb interactions in ionic compounds 23 ) that typically have little angular dependence.A larger radial cutoff is also useful for capturing dispersion interactions in Van der Waals structures 141 .
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 .

Training iterations
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 N P T ensemble.The target pressure is set to zero, and the target temperatures range from 100 K 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: This results in a rapid prediction for the test dataset.The RMSEs for energy, force, and virial are 1.2 meV/atom, 41.6 meV/ Å, and 8.5 meV/atom, respectively.These values are already comparable to those for the training dataset, indicating that we can actually stop here and use NEP-iteration-1 as the final model.However, for added confidence, it is generally advisable to perform at least one more iteration.Therefore, we combine the test dataset (100 structures) with the training dataset (150 structures) to form an expanded training dataset (250 structures), and then train a new model named NEPiteration-2.With this new NEP model, we generate another test dataset with 100 structures, using similar procedure as above but with a simulation time of 10 ns (instead of 10 ps), driven by NEP-iteration-2 for each temperature.The test RMSEs for NEP-iteration-2 are 0.5 meV/atom (energy), 33.5 meV/ Å (force), and 8.9 meV/atom (virial), respectively.Both the energy and force RMSEs are smaller than those for the previous iteration, indicating the improved performance of NEPiteration-2 compared to NEP-iteration-1.
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 NEPiteration-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.

Phonon dispersion relations
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 195 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 196 , 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.

Thermal conductivity from EMD
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: There are three input blocks.In the first block, we specify the NEP potential file and set the initial temperature to 300 K.The second block represents an equilibration run of 500 ps in the N P T ensemble, aiming to reach a target temperature of 300 K and a target pressure of zero.The third block corresponds to a production run of 10 ns in the N V E ensemble, with heat current sampled every 20 steps.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 Wm −1 K −1 from the EMD method.In this work, all statistical errors are calculated as the standard error of the mean.

Thermal conductivity from HNEMD
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: There are also three input blocks, and only the production block differs from the case of EMD.Here, the temperature is controlled using the Nose-Hoover chain thermostat, and an external force in the x direction with F e = 2 × 10 −5 Å−1 is applied.The production run has 10 ns in total.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 54 .
From the HNEMD simulations, we also obtain the spectral thermal conductivity κ(ω) [cf.Eq. (39)].By combining this with the spectral conductance G(ω) [cf.Eq. (37)] in a ballistic NEMD simulation (details provided below), we calculate the phonon MFP spectrum as which is a generalization of Eq. ( 8).The calculated λ(ω) is shown in Fig. 7(c).Remarkably, in the low-frequency limit, λ(ω) can go well beyond one micron.With κ(ω) and λ(ω), one can calculate the spectral apparent thermal conductivity κ(ω, L) according to Eq. ( 10) and obtain the apparent thermal conductivity at any length L using Eq.(11).The results are depicted by the solid line in Fig. 7(d).

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 steadstate 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 sanitychecking the results.

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 197 to be about 250 ± 10 W m −1 K −1 using a Tersoff potential 195 , 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.
According to a recent unpublished study by Wu et al. 181 , the underestimation of thermal conductivity by MLPs could potentially be attributed to small but finite force errors compared to the reference data, leading to extra phonon scatterings.Based on the fact that the force errors form a Gaussian distribution, similar to the random forces in the Langevin thermostat, a method for correcting the force-error-induced underestimation of the thermal conductivity from MLPs is proposed 181 .This correction involves conducting a series of HNEMD simulations with the temperature being controlled by a Langevin thermostat with various relaxation times τ T .Each component of the random force follows a Gaussian distribution with zero mean and a variance of where m is the average atom mass in the system and ∆t is the integration time step.When the random forces in the Langevin thermostat and the force errors in the MLP (with a RMSE of σ mlp at a particular temperature)  178,198,199 and previously (uncorrected and corrected) NEP-MD simulations 181 .
are present simultaneously, a new set of force errors is created, with a larger variance given by according to the properties of Gaussian distribution.After obtaining κ(σ tot ) at different σ tot , the thermal conductivity with zero total force error κ(σ tot = 0) can be obtained from the following relation 181 : where β is a fitting parameter.
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 R bas , N A bas ) 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 181 , which also uses a NEP model [Fig.9(b)].In Fig. 9(b), we also show the results for other temperatures 181 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 181 .
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 MDbased 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.
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 20 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.

FIG. 1 .
FIG. 1. Number of publications (up to March 10th, 2024) on heat transport MD simulations using MLPs as a function of year, with detailed information in Table I and Table II.

FIG. 2 .
FIG. 2. Ballistic-to-diffusive transition of the apparent thermal conductivity κ(L).(a)-(b) a toy model with a single phonon MFP of 1 µm and a diffusive thermal conductivity of κ = 1000 W/mK; (c)-(d) a toy model with two phonon MFPs, one of 0.1 µm, the other 1 µm, with diffusive conductivity of 500 W/mK.The dots in each panel represent a few special lengths, from 0.2 µm to 5 µm.In (a) and (c), the dashed lines represent the ballistic limit.
max and l 5body max as in the 3-body part), see Ref. 25.

FIG. 5 .
FIG. 5. (a) Evolution of RMSEs of energy, force, and virial with respect to training generations (steps).(b) Comparison of force, (c) energy, and (d) virial calculated by NEP against DFT reference data for the initial training dataset.

FIG. 7 .
FIG. 7. Thermal conductivity of crystalline silicon at 300 K from three MD-based methods using the herein developed NEP.(a) Results from 50 independent EMD runs (thin solid lines), along with their average (thick solid line) and error bounds (thick dashed lines); (b) Results from 4 independent HNEMD runs (thin solid lines), along with their average (thick solid line) and error bounds (thick dashed lines); (c) Phonon MFP spectrum calculated using spectral decomposition method; (d) Results from NEMD simulations (red symbols), matching the κ(L) curve from the HNEMD-based formalism.

FIG. 8 .
FIG.8.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. 9. (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 is 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 values178,198,199 and previously (uncorrected and corrected) NEP-MD simulations181 .
Table I and Table II.

TABLE II .
Applications of MLPs in MD simulations of heat transport up to March 10th, 2024.

TABLE III .
32mparison of the energy and force RMSEs and computational speed for MTP, DP (after model compression), and NEP.The computational speed is assessed by running MD simulations for 10 5 steps in the N V T ensemble for a graphene system containing 24 800 atoms, using gpumd (NEP) or lammps32with version 23 Jun 2022 (MTP and DP).For GPU-based tests (DP and NEP), a single Nvidia RTX 3090 is used; for CPU-based tests (MTP), 64 AMD EPYC 7H12 cores are used.