Combining linear-scaling quantum transport and machine-learning molecular dynamics to study thermal and electronic transports in complex materials

We propose an efficient approach for simultaneous prediction of thermal and electronic transport properties in complex materials. Firstly, a highly efficient machine-learned neuroevolution potential (NEP) is trained using reference data from quantum-mechanical density-functional theory calculations. This trained potential is then applied in large-scale molecular dynamics simulations, enabling the generation of realistic structures and accurate characterization of thermal transport properties. In addition, molecular dynamics simulations of atoms and linear-scaling quantum transport calculations of electrons are coupled to account for the electron-phonon scattering and other disorders that affect the charge carriers governing the electronic transport properties. We demonstrate the usefulness of this unified approach by studying electronic transport in pristine graphene and thermoelectric transport properties of a graphene antidot lattice, with a general-purpose NEP developed for carbon systems based on an extensive dataset.


I. INTRODUCTION
Thermal and electronic transports are two fundamental properties of a material.For simple solids, computational methods based on the electron and phonon Boltzmann transport equations [1] have been widely used to compute the transport properties mediated by the heat and charge carriers.There are a handful computational programs available for doing these calculations, such as shengbte [2], phono3py [3], kaldo [4], and gpupbte [5] for thermal transport and epw [6], perturbo [7], and phoebe [8] for electronic transport.However, these methods can only efficiently deal with relatively simple systems and are generally not applicable to complex systems that cannot be properly represented by small periodic supercells.
To efficiently compute transport properties in complex systems one must resort to linear-scaling methods, i.e., methods with the computational cost that scales linearly with respect to the number of atoms in the periodic supercell.For thermal transport, molecular dynamics (MD) simulation is such a linear-scaling method at the atomistic level [9], provided that the interatomic potential used is a classical one and has a finite cutoff.Nowadays, machine-learned potentials (MLPs) [10] have been routinely applied in MD simulations of thermal transport.Particularly, the neuroevolution potential (NEP) [11][12][13] has been developed with a focus on thermal transport applications and has excellent computational efficiency.
For electronic transport, there are also linear-scaling quantum transport (LSQT) methods [14] based on semiempirical tight-binding (TB) models.The electronphonon coupling in LSQT calculations can be captured by the bond-length dependent hopping integrals in the electron TB Hamiltonian [15].This has been explored using either specific phonon dynamics [16,17] or MD simulations [18][19][20].Static-disorder approximation of the electron-phonon coupling has also been used for organic crystals [21,22], graphene [23] and a carbon nanotube [24].Among these, the combined MD-LSQT approach is the most flexible one, but it has not been widely used.The major reason is that there has been no accurate interatomic potential to drive MD simulations for a general system.Another reason is that there is so far no publicly available implementation of this approach.
In this paper, we propose to combine machine-learning molecular dynamics (MLMD), namely, MD driven by a MLP, and LSQT with a bond-length-aware TB model, to study the thermal, electronic, and thermoelectric transport properties of complex materials that are beyond the reach of conventional methods.We call the combined method MLMD-LSQT.For the MLP, we choose to use the highly efficient NEP approach [11][12][13] as implemented in the open-source graphics processing units molecular dynamics (gpumd) package [25].By training against quantum-mechanical density functional theory (DFT) data, a NEP model can be constructed on demand, which can then be used to perform large-scale MD simulations to obtain realistic structures and thermal transport properties.For the LSQT part, we also implement it into the gpumd package (version 3.9) to couple electron and ion motions.To show the usefulness of this unified approach, we construct a generalpurpose NEP for carbon systems and study thermal, electronic, and thermoelectric transport properties of patterned graphene that has large-scale structural features.

II. THE MLMD-LSQT APPROACH
At the core of our method is the NEP approach [11][12][13] for MLP construction.It uses Chebyshev and Legendre polynomials to construct a local atom-environment descriptor of a given atom which is then mapped to the site energy U i of this atom via a feed-forward neural network.The free parameters in the neural network as well as the descriptor are optimized though the minimization of a loss function using an evolutionary algorithm.The loss function is defined as a weighted sum of the root mean square errors (RMSEs) of energy, force, and virial between predictions and DFT target results in combination with regularization terms.This method as implemeted in gpumd [25] has been shown to be able to achieve simultaneously the accuracy of DFT calculations and the computational cost of empirical potentials, allowing for large-scale MD simulations up to 8.1 million atoms using a single 40-gigabyte graphics processing units (GPU) card [26].
The LSQT method can be used to calculate electrical conductivity in large systems, but the prerequisite is to construct an electron Hamiltonian incorporating electron-phonon coupling and other disorders [14].By using a bond-length-aware TB model to configurations generated from MD simulations, electron-phonon coupling and other structural disorders can be effectively described.For dissipative electron transport, there are two equivalent ways to compute the electrical conductivity, one is based on the velocity-auto-correlation and the other is based on the mean-square displacement [14].For the purpose of the present work, we found that the velocity-auto-correlation approach is more convenient because the time intervals used in the calculations are quite small, and the mean-square-displacement approach is only beneficial when the time intervals are large [27].
In the velocity-auto-correlation approach, the electrical conductivity at energy E and correlation time t can be calculated as an integral where e is the elementary charge, Ω is the system volume, Ĥ is the electron Hamiltonian operator, δ(E − Ĥ) is the energy resolution operator, V is the velocity operator, and V (τ ) = e i Ĥτ V e −i Ĥτ is the time-evolved velocity operator.To facilitate the discussion, we denote the trace in the integral as C(E, τ ).The coupled MLMD-LSQT algorithm can be represented as follows: 1. Starting from an initial structure, run MLMD for a number of steps in the isothermal or isothermalisobaric ensemble to achieve equilibrium.

Perform MLMD simulation for a number of steps:
(a) Evolve the atomic system from step n − 1, {r i (n∆t − ∆t)}, to step n, {r i (n∆t)}, by a time step of ∆t according the NEP interatomic potential.
(b) Calculate the electron Hamiltonian and velocity operators at step n according to the atom positions {r i (n∆t)}.
(c) Calculate C(E, n∆t) using the electron Hamiltonian at the current step.In this step, linear-scaling techniques [14], including sparse matrix-vector multiplication, random phase approximation of trace, Chebyshev expansion of quantum evolution operator, and kernel polynomial method [28] for energy resolution operator, are used.
After obtaining C(E, τ ) at a number of discrete time points, it can be numerically integrated to calculate the electrical conductivity according to Eq. 1.This approach was implemented into the gpumd package and was available starting from version 3.9.Besides, the electronic density of states (DOS) was also implemented according to the following expression:

III. CASE STUDY OF A GRAPHENE ANTIDOT LATTICE
As a proof of concept, we apply the MLMD-LSQT approach to study the thermoelectric transport in a GAL [29], also known as graphene nanomesh [30], a graphene sheet with patterned holes.Thermoelectric effects in graphene nanostructures have been extensively studied, and GALs have been identified as one of the promising candidates for good thermoelectric materials [31].However, previous works have only studied the ballistic electronic transport regime [32][33][34], without considering finite-temperature effects.
Fig. 1 shows the atomistic structure of the system under investigation.The simulation domain cell of the GAL sample contains 187 200 atoms and has a dimension of about 88.5 nm × 76.7 nm in the xy-plane, which can be considered as a two-dimensional (2D) system when periodic boundary conditions are applied to the in-plane directions.The thickness of the system was taken as 0.335 nm in calculating the volume.The primitive cell for the GAL contains 156 atoms, a complexity that challenges conventional numerical methods based on the electron and phonon Boltzmann transport equations.However, this kind of complex structures are well-suited for the MLMD-LSQT approach.To construct the Hamiltonian and velocity operators, we employed a p z -orbital TB model with a bond-length dependent hopping parameter where t 0 = −2.7 eV, r 0 = 1.42 Å, and r ij is the distance between the atom pair i and j.The model with a fixed hopping parameter t 0 has been used in previous works [35,36] that did not account for electron-phonon coupling.The real-space Hamiltonian and velocity (assuming to be in the x direction) operators can be written as where x i is the x-position of atom i.

A. Training a general-purpose NEP for carbon systems
Although for the scope of the current work, it suffices to train a specialized NEP model for GAL, it is our broader objective to train a general-purpose carbon potential based on the extensive dataset as used for constructing a Gaussian approximation potential [37].Using this dataset and the hyperparameters given in Appendix A, we trained a general-purpose NEP model for carbon The seemingly large RMSE values are typical for generalpurpose carbon systems, as similar ones were reported in or can be extracted from previous works [13,[37][38][39][40].

B. Thermal transport
For a complete study of thermoelectric transport, the lattice (phonon) thermal conductivity κ ph must be evaluated.To this end, we calculated κ ph for the GAL model with 187 200 atoms using the HNEMD method [41].In this method, an external driving force is exerted on each atom i, driving the system out of equilibrium.Here, F e is the driving force parameter with the dimension of inverse length and r ij ≡ r j −r i , r i being the position of atom i.After a steady state is achieved, the lattice thermal conductivity tensor κ αβ ph can be computed from the relation where T is the system temperature, Ω is the system volume, and ⟨J α ⟩ is the ensemble average of the heat current [42] In this case study, we only consider the condition of 300 K and zero in-plane pressure.The input script for gpumd is given in Appendix B. The time convergence of κ ph is shown in Fig. 3(a).We have performed 4 independent simulations in both the x and the y directions and averaged the results over the two directions as the system is essentially isotropic.For the GAL in our case study, the simulation temperature (300 K) is much lower than the Debye temperature (on the order of 2000 K), and the classical MD simulation thus significantly overestimates the modal heat capacity, which in turn leads to an overestimation of the thermal conductivity.Fortunately, there exists a feasible correction for the missing quantum statistics, as has been successfully applied to amorphous [43,44] and fluid [45] systems described by NEP models.In this quantum correction method, the spectral thermal conductivity κ(ω) as calculated from the HNEMD method [41] is multiplied by a quantum-to-classical factor p(x) = x 2 e x / (e x − 1) 2 , where x = ℏω/k B T , ω is the phonon frequency, ℏ is the reduced Planck constant, and k B is the Boltzmann constant.The spectral decomposition of κ ph corresponding to the classical results is depicted in Fig. 3(b), where the quantum-corrected results are also shown.The classical value of κ ph is 13.3 W/(m K), which becomes 6.7 W/(m K) after quantum correction.The quantum corrected value will be used later.

C. Electronic and thermoeletric transports
The time-dependent electrical conductivity Σ(E, t) (see Appendix C for the input script and more calculation details) converges in the diffusive transport regime, and one can obtain the so-called semi-classical electrical conductivity Σ(E) by averaging Σ(E, t) over a proper range of correlation time, According to Fig. 4(a), it is a good choice to set t 1 = 80 fs and t 2 = 100 fs.The semi-classical electrical conductivity Σ(E) can be regarded as the transport distribution function (TDF) [46][47][48][49][50][51] for thermoelectric transport.The calculated TDF as well as DOS are presented in Fig. 4(b).The anti-dots induce a considerable band gap of about 0.8 eV.This band gap is then also the transport gap.
From the TDF, we then calculated the transport coefficients at 300 K for a range of chemical potential µ.We first define the functionals (n = 0, 1, 2) of the TDF: where is the Fermi-Dirac distribution.The electrical conductivity σ(µ, T ), Seebeck coefficient S(µ, T ), and electronic thermal conductivity κ el (µ, T ) can be expressed in terms of these functionals as σ(µ, T ) = X 0 (µ, T ); ( 12) The calculated results are presented in Fig. 4(c)-(e).The finite-temperature electrical conductivity σ(µ, T ) resembles the TDF, but with smearing resulting from the Fermi-Dirac distribution.The Seebeck coefficient has a negative peak for electrons and a positive peak for holes.The electronic thermal conductivity resembles the electrical conductivity in shape which is in line with the Wiedemann-Franz law.
Based on these transport coefficients and the (quantum-corrected) phonon thermal conductivity κ ph , one can define the dimensionless figure of merit as Due to the competition between the various transport coefficients, zT develops peaks for both electron and hole transport, at µ = 0.47 eV and −0.44 eV, respectively.The transport is asymmetric between electron and hole, showing a maximum zT = 0.47 for hole and a maximum zT = 0.16 for electron.Experimentally, thermoelectric transport properties have been measured for single-and bi-layer graphene nanomeshes with the neck width down to 8 nm [52].This neck width between the nearest antidot pairs is a few time larger than that we studied and the measured thermal conductivity values (of the order of 100 W/(m K)) are significantly larger than our prediction.On the other hand, there are geometrical disorders in the experimental samples, namely, variations in the positions and sizes of the antidots, which, according to previous calculations [35,36], can lead to suppressed electrical conductivity.Therefore, the relatively high zT values we predicted remain a challenge for experimental realization.

IV. SUMMARY AND CONCLUSIONS
In summary, we have introduced a numerical approach for simultaneous prediction of thermal and electronic transport properties in complex materials.This approach, based on MLMD and LSQT, offers an excellent efficiency with a computational cost that scales linearly with the system size.For a given material, a highly efficient NEP is first constructed on demand.This MLP can be used to perform large-scale MD to obtain realistic structures and accurate thermal transport properties.By combining the time-evolution of electrons and atoms during the MD simulation, electron-phonon scattering and other disorders for the charge carriers can be naturally captured and the various electronic transport properties can be obtained.
As an illustrative example, we have investigated the thermoelectric transport properties of a type of graphene antidot lattices (GALs), predicting its relatively high thermoelectric efficiency at room temperature.We recognize the necessity of future work to conduct a more comprehensive study of the thermoelectric transport in GALs using our proposed approach.Subsequent research endeavors may consider integrating machine-learning techniques to explore the vast design space, similar to the methods employed in thermal transport studies in these systems [53,54].• The first parameter x means that the transport is along the x direction.We have calculated 10 times along the x directions and also 10 times along the y direction and averaged the results.
• The second parameter refers to the number of Chebyshev moments in the kernel-polynomial method [28] for both the DOS and conductivity calculations.A value of 3000 is large enough here.
• The next three parameters are respectively the number of energy points to be considered, the minimum energy and the maximum energy.Here we calculated the transport properties from −8.1 eV to 8.1 eV, with an interval of 1.62 meV.
• The last parameter is an energy threshold that needs to be larger than the energy range of the tight-binding model.Here, a value of 8.2 eV is sufficient.This parameter can be determined by a trial-and-error approach. #

FIG. 1 .
FIG. 1.(a) Atomistic structure of an example graphene antidot lattice (GAL) system studied in this work.(b) Illustration of the primitive cell containing 156 atoms, enclosed by the parallelogram.
systems.The training results are shown in Fig. 2.After a few hundred thousand training steps, the RMSEs of energy, force, and virial all converge [Fig.2(a)], and their converged values are 45 meV/atom, 599 meV/ Å, and 105 meV/atom, respectively.The predicted data are compared to the DFT reference ones in Fig. 2(b)-(d).

FIG. 2 .FIG. 3 .
FIG. 2. (a) Evolution of the energy, force, and virial loss values as a function of the number of training generations for a general-purpose NEP for carbon systems based on an extensive dataset [37].(b)-(d) Comparison between NEP predictions and DFT reference values for energy, force, and virial.The RMSE values are indicated in each panel.

ACKNOWLEDGMENTSZ.
Fan and H. Dong were supported by the National Natural Science Foundation of China (NSFC) (No. 11974059) and the Research Fund of Bohai University (No. 0522xn076).P. Ying was supported by the Israel Academy of Sciences and Humanities & Council for Higher Education Excellence Fellowship Program for International Postdoctoral Researchers.