Accurate prediction of heat conductivity of water by a neuroevolution potential

TL;DR

This paper introduces a machine learning-based potential combined with advanced simulation techniques to accurately predict water's heat conductivity, aligning well with experimental data across various conditions.

Contribution

It develops a neuroevolution-based potential and integrates Green-Kubo and spectral methods to incorporate quantum effects in heat conductivity predictions.

Findings

Achieves quantum-mechanical accuracy in heat conductivity predictions.

Matches experimental results across temperature and pressure conditions.

Provides a computationally efficient alternative to quantum simulations.

Abstract

We propose an approach that can accurately predict the heat conductivity of liquid water. On the one hand, we develop an accurate machine-learned potential based on the neuroevolution-potential approach that can achieve quantum-mechanical accuracy at the cost of empirical force fields. On the other hand, we combine the Green-Kubo method and the spectral decomposition method within the homogeneous nonequilibrium molecular dynamics framework to account for the quantum-statistical effects of high-frequency vibrations. Excellent agreement with experiments under both isobaric and isochoric conditions within a wide range of temperatures is achieved using our approach.

Tables

Table 1. Table 1: Calculated classical and quantum-corrected heat conductivity values (in units of Wm -1 K -1 ) at different temperatures T 𝑇 T (in units of K) under both isobaric and isochoric conditions. The numbers within the parentheses are statistical uncertainties for the last significant number(s).

	isobaric		isochoric
$T$	classical	quantum	classical	quantum
275	0.91(4)	0.58(5)	0.89(6)	0.57(10)
287.5	0.88(3)	0.59(4)	0.89(4)	0.58(5)
300	0.89(5)	0.60(7)	0.89(5)	0.60(7)
312.5	0.90(3)	0.62(4)	0.82(5)	0.58(7)
325	0.88(5)	0.64(6)	0.94(8)	0.66(11)
337.5	0.85(6)	0.62(7)	0.90(4)	0.67(5)
350	0.84(5)	0.66(6)	0.89(6)	0.69(7)
362.5	0.86(5)	0.68(6)	0.93(4)	0.71(5)
375	0.83(8)	0.66(9)	0.90(6)	0.72(8)
387.5	0.84(6)	0.67(6)	0.88(7)	0.75(9)
400	0.84(6)	0.70(7)	0.90(8)	0.74(9)
412.5	0.81(6)	0.66(7)	0.89(4)	0.76(4)
425	0.78(4)	0.64(5)	0.89(4)	0.76(5)
437.5	0.76(5)	0.64(6)	0.96(5)	0.81(6)
450	0.77(5)	0.66(7)	0.95(5)	0.80(6)
462.5	0.74(5)	0.64(6)	0.91(4)	0.80(5)
475	0.75(5)	0.64(6)	0.96(6)	0.85(7)
487.5	0.72(7)	0.61(8)	0.90(9)	0.80(12)
500	0.72(5)	0.64(6)	0.99(6)	0.85(7)

Equations

$${\kappa }^{\mathrm{total}}(t)=\frac{1}{3k_{\mathrm{B}}V{T}^2}{\int }_0^{t}d{t}^{\prime }⟨\mathbf{J}(t^{\prime })\cdot \mathbf{J}(0)⟩,$$

$$\mathbf{J}={\mathbf{J}}^{\mathrm{k}}+{\mathbf{J}}^{\mathrm{p}}.$$

$${\mathbf{J}}^{\mathrm{k}}=\sum _i{\mathbf{v}}_i{E}_i,$$

$${\mathbf{J}}^{\mathrm{p}}=\sum _i\sum _j{\mathbf{r}}_{ij}\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}\cdot {\mathbf{v}}_i.$$

$${\kappa }^{\mathrm{total}}(t)={\kappa }^{\mathrm{pp}}(t)+{\kappa }^{\mathrm{kk}}(t)+{\kappa }^{\mathrm{pk}}(t),$$

$$p(x)=\frac{x^2{e}^x}{(e^x-1{)}^2},$$

$${\mathbf{F}}_i={\mathbf{F}}_{\mathrm{e}}\cdot \sum _{j\ne i}{\mathbf{r}}_{ij}\otimes \frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}$$

$${\kappa }^{\mathrm{pp}}(\omega )=\frac{2}{VT{F}_{\mathrm{e}}}{\int }_{-\infty }^{\infty }\text{d}t{e}^{\text{i}\omega t}\sum _i\sum _{j\ne i}⟨x_{ij}\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}(0)\cdot {\mathbf{v}}_i(t)⟩.$$

$${\kappa }^{\mathrm{pp}}={\int }_0^{\infty }\frac{\text{d}\omega }{2\pi }{\kappa }^{\mathrm{pp}}(\omega ).$$

$${\kappa }_{\mathrm{q}}^{\mathrm{pp}}(\omega )={\kappa }^{\mathrm{pp}}(\omega )p(x).$$

$${\kappa }_{\mathrm{q}}^{\mathrm{pp}}={\int }_0^{\infty }\frac{\text{d}\omega }{2\pi }{\kappa }_{\mathrm{q}}^{\mathrm{pp}}(\omega ).$$

$${\kappa }_{\mathrm{q}}^{\mathrm{total}}={\kappa }_{\mathrm{q}}^{\mathrm{pp}}+{\kappa }^{\mathrm{kk}}+{\kappa }^{\mathrm{pk}}.$$

$$\kappa (\omega )=c(\omega )v^2(\omega )\tau (\omega ),$$

Full text

I Introduction

Heat transport in fluids involves both interatomic interactions and diffusion, and classical molecular dynamics (MD) simulation is a viable method for computing heat conductivity by including both the interaction and diffusion contributions. Extensive MD simulations [1, 2, 3, 4, 5, 6, 7] have been performed to calculate the heat conductivity of water using various empirical force fields, such as SPC/E [8], TIP4P [9], TIP4P/2005 [10], and ReaxFF [11]. However, The force fields were found to have a major influence on the calculated thermal conductivity and quantitative agreement between simulations and experimental measurements in a wide range of temperatures has not been achieved for any force field so far.

One of the reasons for the disagreement between computations and measurements is the inaccuracy of the empirical force fields. Although classical MD simulations of heat transport can also be driven by interactions computed by quantum-mechanical density functional theory (DFT) [12, 13], this approach is currently not efficient enough and has not been extensively applied to heat transport in water. Recently, machine-learned potentials (MLPs) [14] have emerged as an alternative that can achieve the accuracy of DFT with a small fraction of the cost. Recent studies [15, 16, 17, 18, 19] have demonstrated the high accuracy of MLPs in modeling the thermodynamics of water in various phases. The linear-scaling computational cost with respect to the number of atoms enabled efficient MD simulation of heat transport in complex systems that is beyond the reach of perturbative methods [20, 21, 22, 23, 24, 25]. A deep potential (DP) model [26] has been developed to calculate the heat conductivity of water in a range of temperatures. However, only a qualitative agreement with experiments has been achieved. It has been not clear if MLPs can reliably predict the heat conductivity of liquid water at a wide range of thermodynamic conditions.

In this work, we developed a MLP for water within the neuroevolution potential (NEP) framework [27, 28, 29], which is an efficient MLP framework that has been developed with a particular emphasis on heat transport applications. The accuracy of the developed NEP model is demonstrated by the radial and angular distribution functions as compared to DFT results. We performed equilibrium molecular dynamics (EMD) simulations to calculate the heat conductivity using the Green-Kubo relation [30, 31]. The results from classical MD simulations driven by NEP do not match experiments quantitatively. However, by applying a quantum-statistical correction based on the spectral heat conductivity computed within the homogeneous nonequilibrium molecular dynamics (HNEMD) approach [32], a quantitative agreement with experiments can be achieved for a wide range of temperatures at both isobaric and isochoric conditions.

II A NEP model for liquid water

The NEP approach as implemented in the gpumd package [33] has been introduced in Ref. 27 and improved later [28, 29]. This approach follows the work of Behler and Parrinello [14] to model the site energy of an atom as an artificial neural network (ANN), where the input layer consists of a descriptor vector of high dimensions. The descriptor components are invariant with respect to the translation, rotation, and permutation of atoms of the same kind. For explicit expressions of the descriptor components, we refer to Ref. 29. The name NEP comes from the algorithm for training the ANN, which is a separable natural evolution strategy (SNES) [34].

To train a NEP model for liquid water, we used the “refinement” data set for liquid water taken from Ref. 19, computed at the quantum-mechanical DFT level with the strongly constrained and appropriately normed (SCAN) functional [35]. There are 1888 structures (each with 128 H2O molecules) in total, and we randomly selected 1388 for training and 500 for testing. For more details on the generation of the reference data, we refer to Ref. 19.

The hyperparameters we used in the NEP model are as follows. The NEP descriptor consists of a number of radial and angular components [27, 28, 29]. For the radial components, we used a cutoff radius of 6 Å and ten radial functions (each being a linear combination of 10 basis functions). For the angular components, we used a cutoff radius of 4 Å, eight radial functions (each being a linear combination of 8 basis functions), three-body correlations up to $l=4$ in the spherical harmonics, and four-body correlations up to $l=2$. The ANN in the NEP model has a single hidden layer and we used 100 neurons for this layer.

We trained the NEP model for 300,000 generations using the SNES algorithm and the loss terms for energy, force, and virial in the test set are largely converged [Fig. 1(a)]. The predicted energy, force, and virial for the test set are compared to the DFT reference data in Figs. 1(b) to 1(d), showing good correlations. Quantitatively, the root mean square errors (RMSEs) for energy, force, and virial are 0.89 meV/atom, 76 meV/Å, and 5.2 meV/atom in the training data set, and are 1.0 meV/atom, 73 meV/Å, and 5.0 meV/atom in the test data set. The level of accuracy is comparable to those reported in previous works on identical or similar data sets [19, 26].

Our NEP model can achieve not only a high accuracy but also a high computational speed. To show this, we compare the computational speeds of NEP as implemented in the gpumd package (version 3.6) [33] and the SPC/E force field as implemented in the lammps package (the 23 Jun 2022 version) [36]. For the SPC/E force field, the Coulomb interactions were evaluated using the particle-particle particle-mesh (PPPM) method with a real-space cutoff distance of 12 Å and a relative accuracy of $10^{-6}$ in force calculations. Figure 2 shows that our NEP model is literally as fast as the SPC/E force field for comparable amounts of computational resources.

To validate the accuracy of the trained NEP model in MD simulations, we compare the radial distribution function (RDF) and angular distribution function (ADF) obtained by classical MD simulations driven by the NEP model and DFT calculations, both at 300 K and 1 bar. As shown in Fig. 3, good agreement is achieved for the RDFs for O-O pairs, $g_{\rm{OO}}(r)$, and O-H pairs, $g_{\rm{OH}}(r)$, and the ADF for O-O-O triplets $g_{\rm{OOO}}(\theta)$. The DFT results were obtained using a small cell with 384 atoms, while the NEP results were obtained using a much larger cell with 3000 atoms, which explains the much smoother distribution functions from the NEP model. The good agreement here indicates that our NEP model can accurately reproduce the dynamics of liquid water, which is a prerequisite for the reliable study of heat transport.

III Heat conductivity of liquid water from MD simulations

III.1 Classical heat conductivity of liquid water

We calculated the heat conductivity of liquid water using the well established Green-Kubo method [30, 31], in which the running heat conductivity $\kappa^{\rm total}(t)$ (the meaining of the superscript “total” will become clear below) is calculated as a time integral of the heat current autocorrelation function (HCACF) $\langle\bm{J}(t)\cdot\bm{J}(0)\rangle$:

[TABLE]

where $k_{\rm B}$ is Boltzmann’s constant, $T$ and $V$ are the temperature and volume of the system, respectively. The heat current $\bm{J}(t)$ is sampled at an equilibrium state. For liquid molecules, the heat current has two contributions:

[TABLE]

The kinetic term (also called convective term) is

[TABLE]

and the potential term for many-body potentials such as our NEP model is [37]

[TABLE]

Here, $E_{i}=\frac{1}{2}m_{i}\bm{v}_{i}^{2}+U_{i}$ is the total energy of atom $i$, where $m_{i}$, $\bm{v}_{i}$, and $U_{i}$ are respectively the mass, velocity, and potential energy of atom $i$. According to the decomposition of the heat current, the heat conductivity can be decomposed into three terms:

[TABLE]

where the potential-potential term $\kappa^{\rm pp}$, the kinetic-kinetic term $\kappa^{\rm kk}$, and the cross term $\kappa^{\rm pk}$ correspond to the following HCACFs: $\langle\bm{J}^{\rm p}(t)\cdot\bm{J}^{\rm p}(0)\rangle$, $\langle\bm{J}^{\rm k}(t)\cdot\bm{J}^{\rm k}(0)\rangle$, and $\langle\bm{J}^{\rm p}(t)\cdot\bm{J}^{\rm k}(0)\rangle+\langle\bm{J}^{\rm k}(t)\cdot\bm{J}^{\rm p}(0)\rangle$.

The Green-Kubo method is based on EMD, where the system is first equilibrated in the $NVT$ (constant number of atoms $N$, constant volume $V$, and constant target temperature $T$) or $NpT$ (constant target pressure $p$) ensemble to reach an equilibrium state and the heat currents are then sampled in the $NVE$ (constant energy $E$) ensemble. In all the MD simulations in this work, we used a time step of 0.1 fs, which has been tested to be small enough. In the EMD simulations, we used an equilibration time of 50 ps and a production time of 10 ps. For each thermodynamic state with a given temperature and pressure (or density), we performed about 50 independent runs and calculated the statistical error as the standard error between the independent runs. We have tested the effects of finite simulation cells and found that the heat conductivity is essentially unchanged when the linear size of a cubic cell increases from 3 to 11 nm. We chose to use a cell with a linear size of about 6 nm containing 24576 atoms (8192 water molecules) for all the subsequent calculations.

Figures 4(a)-4(c) show the running heat conductivity components $\kappa^{\rm pp}(t)$, $\kappa^{\rm kk}(t)$, and $\kappa^{\rm pk}(t)$, respectively. In the interval from $t=3$ to 5 ps, all the components show stable oscillations only, without an overall increasing or decreasing trend. We therefore average the running heat conductivity over this time interval for each independent run. With about 50 independent runs, we thus obtained a mean value of each heat conductivity component and a statistical error estimate. These are shown in Fig. 4(d). Among the three components, $\kappa^{\rm kk}$ is about one order of magnitude smaller than $\kappa^{\rm pp}$ and $\kappa^{\rm pk}$ is essentially zero.

Using the Green-Kubo method, we computed the total heat conductivity $\kappa^{\rm total}$ of liquid water from 275 to 500 K (with a step of 12.5 K) under both isobaric and isochoric conditions. Isobaric conditions were achieved by using the $NpT$ ensemble [38] with a target pressure of 30 bar. Isochoric conditions were achieved by using the $NVT$ ensemble [39] with a fixed density of 1 g/cm3. Our results for isobaric and isochoric conditions are presented in Figs. 5 and  6, respectively, along with the experimental data from national institute of standards and technology (NIST) [40, 41] and previous theoretical ones. In the isobaric case, our heat conductivity values (the classical ones) are very close to the previous ones obtained by using the DP approach [26], except for a noticeable difference around 300 K. However, the heat conductivity values from both DP and our NEP are significantly higher than the experimental ones, particularly at the lower temperatures. The predicted results from the empirical force fields (SPC/E, TIP4P, and TIP4P/2005) [6] show a more complex pattern: they are relatively high at about 400 K but can be close to or lower than experimental values at both the low- and high-temperature limits. Because the experimental data from NIST [40, 41] have small uncertainties (a few percent at most), the results in Fig. 5 show that all the theoretical predictions do not quantitatively agree with the experiments. As we will argue below, nuclear quantum effects (NQEs) play an important role here.

III.2 Quantum-corrected heat conductivity of liquid water

In classical MD simulations, the vibrations in the system follow the classical statistics, with all degrees of freedom being fully activated regardless of the temperature and frequency. However, according to quantum statistics, high-frequency degrees of freedom are frozen at low temperatures. Quantitatively, a degree of freedom with frequency $\omega$ at temperature $T$ is only activated with the following probability:

[TABLE]

where $x=\hbar\omega/k_{\rm{B}}T$, $\hbar$ being the reduced Planck constant. Frequency domain quantum correction [42] based on the vibrational density of states (VDOS) has been successfully applied to correct thermodynamic quantities (such as heat capacity) calculated using classical MD. To our best knowledge, this type of quantum correction has not been applied to heat transport in liquid water. Similar to the quantum correction of heat capacity in water [42] based on spectral analysis, heat conductivity can be quantum-corrected based on a spectral heat conductivity, as has been recently demonstrated for amorphous silicon in the context of MLP [25]. Such spectral heat conductivity can be conveniently obtained in the framework of the HNEMD method as developed in Ref. 32.

In the HNEMD simulations, a driving force

[TABLE]

was applied to each atom $i$ of the system to drive the system into a nonequilibrium steady state, in which the classical spectral heat conductivity is calculated. The vector $\bm{F}_{\rm e}$ represents the driving force parameter that is of the dimension of inverse length. For more details on the HNEMD method for many-body potentials, we refer to Ref. 32. Five independent runs, each with a production time of 100 ps, were performed to calculate the mean value and statistical error of the heat conductivity. As only the potential-potential part of the HCACF involves high frequencies that require a quantum correction, here we only apply the HNEMD method to calculate the spectral heat conductivity $\kappa^{\rm pp}(\omega)$, which can be expressed as [32]:

[TABLE]

Here we have assumed that heat transport is along the $x$ direction. The magnitude $F_{\rm e}$ of the driving force parameter is set to $0.001$ Å*-1*, which is sufficiently small to keep the system within the linear-response regime.

Eq. (8) represents a Fourier transform in which the integral is formally from $-\infty$ to $\infty$. In numerical calculations, it is evaluated based on discrete Fourier transform (or more exactly, discrete cosine transform). In MD simulation, the virial-velocity time correlation function $\langle\cdots\rangle$ in Eq. (8) is evaluated at discrete times and only up to a finite upper limit $t_{\rm max}$. A Hann window function is applied before performing the discrete cosine transform. According to Nyquist sampling theorem, $t_{\rm max}$ determines the frequency resolution that can be achieved: a larger value of $t_{\rm max}$ results in a finer frequency resolution. In our calculations, we used $t_{\rm max}=250$ fs, which gives a frequency resolution of $1/2t_{\rm max}$ = 2 THz. This is sufficient for our purpose. Using larger $t_{\rm max}$ does not affect any of our results significantly.

The integration of the spectral heat conductivity over the frequency $\omega$ is $\kappa^{\rm pp}$:

[TABLE]

For the potential-potential part, HNEMD and EMD give consistent heat conductivity, as shown in Fig. 4(d). With the classical spectral heat conductivity available, we can then obtain a quantum-corrected spectral heat conductivity $\kappa_{\rm q}^{\rm pp}(\omega)$ by multiplying $\kappa^{\rm pp}(\omega)$ with the probability $p(x)$:

[TABLE]

Figure 7(a) shows that the quantum correction is significant at 300 K.

After applying this quantum correction, we can obtain the overall quantum-corrected potential-potential part of the heat conductivity as

[TABLE]

Because the kinetic-kinetic and potential-kinetic parts do not involve high-frequency vibrations, they do not need to be quantum corrected. Therefore, we can obtain the quantum-corrected total heat conductivity as

[TABLE]

The total heat conductivity values before and after the quantum correction are listed in Table 1 and are also shown in Figs. 5 and 6.

The quantum-corrected total heat conductivity values agree excellently with experiments in the whole temperature range for both isobaric (Fig. 5) and isochoric (Fig. 6) conditions. Quantitative agreement with experiments for the different temperature and pressure (density) conditions cannot be an accident and it strongly suggests the reliability of our NEP model and the effectiveness of the quantum correction based on the spectral heat conductivity. Particularly, our approach has correctly predicted the much larger heat conductivity under the isochoric condition than that under the isobaric condition at 500 K. This is intuitively understandable as the density for the isochoric condition is higher than that for the isobaric condition, resulting in stronger interatomic interactions that can enhance the potential-potential part of the heat conductivity. Figure 7(b) further shows that this enhancement mainly comes from the vibrations with $\omega/2\pi<10$ THz.

We can now better interpret the theoretical predictions from the DP [26] model and the three empirical force fields [6] as shown in Fig. 5. If the quantum correction were also applied to these predictions, we expect that the DP results would agree well with experiments as well, except for the temperatures close to 300 K. On the other hand, all the three empirical force fields would significantly underestimate the experimental results at 300 K.

IV Discussion

Good agreement between our theoretical calculations and experiments clearly depends on both the accuracy of the NEP model and the effectiveness of the HNEMD-based quantum correction method. Here we discuss the rationales of the HNEMD-based quantum correction method and related approaches.

We start our discussion by presenting a general expression of the spectral heat conductivity:

[TABLE]

where $c(\omega)$ is the modal heat capacity, $v(\omega)$ is the modal group velocity, and $\tau(\omega)$ is the relaxation time of the heat carriers, which are not necessarily phonons but are related to the collective vibrations in the system.

Clearly, one of the NQEs is related to the modal heat capacity: classical MD overestimates the modal heat capacity by exciting any vibrational mode regardless of its frequency and temperature. The quantum correction for this is simple, which is to multiply $c(\omega)$ by $p(x)$ as defined in Eq. (6).

There are essentially no NQEs in the group velocity but there can be complicated NQEs in the relaxation time. This is the case for crystals as has been discussed in the context of phonon Boltzmann transport equation [43, 44, 45]. For example, if classical MD overestimates the population of a given phonon frequency, it leads to overestimated scattering to other phonons. Whether classical MD leads to overestimated or underestimated heat conductivity thus depends on the competition between the NQEs on $c(\omega)$ and $\tau(\omega)$. Even though classical MD might lead to the correct total heat conductivity, the spectral heat conductivity can significantly deviate from the quantum result. Therefore, there is so far no feasible quantum-correction method for heat conductivity of crystals for which phonon-phonon scattering is the major source of resistivity [46]. Particularly, the temperature-rescaling method [47, 48, 49] based on equating the classical and quantum energies has been shown to be infeasible [43, 45].

The situation is different for disordered materials, where the elastic scattering for the vibrational modes by disorder dominates and the population of vibrations has negligible effects on the elastic scattering processes. Therefore, the major quantum effects in classical simulation of disordered systems are from the overestimated modal heat capacity. In this case, the spectral heat conductivity can be quantum corrected by multiplying it with $p(x)$, which is consistent with the HNEMD-based quantum-correction method. The effectiveness of this quantum-correction method has been recently demonstrated for amorphous materials [25] and our current work extends its applicability to liquids by complementing it with EMD simulations for convective heat transport.

While we have only studied liquid water in this work, we believe that our approach is also applicable to other fluids with light elements at relatively low temperatures. However, a more systematic study is needed to evaluate the effectiveness of our approach in other systems. We note that quantum MD methods such as linearized semiclassical initial value representation, centroid MD, and ring-polymer MD have been used to study heat transport of both liquids [50, 51, 52] and solids [53] to account for the NQEs. The relative performance of our approach compared to these quantum MD methods in predicting heat conductivity remains to be explored. Particularly, our approach does not account for NQEs in $\kappa^{\rm kk}$, which is essentially a zero-frequency property similar to the diffusion coefficient. As it has been shown that there are large NQEs in the diffusion coefficient of liquid water [54], we expect that there are also NQEs in $\kappa^{\rm kk}$. However, we note that the classical value of $\kappa^{\rm kk}$ only contributes about 10% to the total heat conductivity. Therefore, even though we have not quantum-corrected $\kappa^{\rm kk}$, we have only ignored little NQEs for the total heat conductivity.

V Conclusions

In summary, we have constructed a NEP model for liquid water that can accurately reproduce structural properties as determined by quantum-mechanical DFT calculations. The NEP model is as efficient as empirical force fields of liquid water in large-scale MD simulations. Heat conductivity values calculated using the Green-Kubo method within classical MD simulations were found to be overestimated against experimental results, particularly for relatively low temperatures. This led us to identify the importance of NQEs in determining the heat conductivity of water. We then proposed a scheme of quantum correction based on the spectral heat conductivity as calculated within the framework of HNEMD simulations, which leads to excellent agreement with experiments under both isobaric and isochoric conditions within a large range of temperatures.

Acknowledgements.

K.X, Y.H, and J.W acknowledge support from the National Natural Science Foundation of China (NSFC) (No. 12172314, 11772278, and 11904300), the Jiangxi Provincial Outstanding Young Talents Program (No. 20192BCBL23029), the Fundamental Research Funds for the Central Universities (Xiamen University: No. 20720210025), and the 111 project (B16029). Z.F. acknowledges support from NSFC (No. 11974059).

Conflict of Interest

The authors have no conflicts to disclose.

Data availability

The source code and documentation for gpumd are available at https://github.com/brucefan1983/GPUMD and https://gpumd.org, respectively. The training and testing results for the NEP model are freely available at https://gitlab.com/brucefan1983/nep-data.