Search

Path-integral molecular dynamics (PIMD) simulations are crucial for accurately capturing nuclear quantum effects in materials. However, their computational intensity and reliance on multiple software packages often limit their applicability at large scales. Here, we present an integration of PIMD methods, including thermostatted ring-polymer molecular dynamics (TRPMD), into the open-source GPUMD package, combined with highly accurate and efficient machine-learned neuroevolution potential (NEP) models. This approach achieves almost the accuracy of first-principles calculations with the computational efficiency of empirical potentials, enabling large-scale atomistic simulations that incorporate nuclear quantum effects. We demonstrate the efficacy of the combined NEP-PIMD approach by examining various thermal properties of diverse materials, including lithium hydride (LiH), three porous metal-organic frameworks (MOFs), and elemental aluminum. For LiH, our NEP-PIMD simulations successfully capture the isotope effect, reproducing the experimentally observed dependence of the lattice parameter on the reduced mass. For MOFs, our results reveal that achieving good agreement with experimental data requires consideration of both nuclear quantum effects and dispersive interactions. For aluminum, the TRPMD method effectively captures thermal expansion and phonon properties, aligning well with quantum mechanical predictions. This efficient NEP-PIMD approach opens new avenues for exploring complex material properties influenced by nuclear quantum effects, with potential applications across a broad range of materials., Comment: 14 pages, 6 figures in the main text; 1 table and 5 figures in the SI

Disordered forms of carbon are an important class of materials for applications such as thermal management. However, a comprehensive theoretical understanding of the structural dependence of thermal transport and the underlying microscopic mechanisms is lacking. Here we study the structure-dependent thermal conductivity of disordered carbon by employing molecular dynamics (MD) simulations driven by a machine-learned interatomic potential based on the efficient neuroevolution potential approach. Using large-scale MD simulations, we generate realistic nanoporous carbon (NP-C) samples with density varying from $0.3$ to $1.5$ g cm$^{-3}$ dominated by sp$^2$ motifs, and amorphous carbon (a-C) samples with density varying from $1.5$ to $3.5$ g cm$^{-3}$ exhibiting mixed sp$^2$ and sp$^3$ motifs. Structural properties including short- and medium-range order are characterized by atomic coordination, pair correlation function, angular distribution function and structure factor. Using the homogeneous nonequilibrium MD method and the associated quantum-statistical correction scheme, we predict a linear and a superlinear density dependence of thermal conductivity for NP-C and a-C, respectively, in good agreement with relevant experiments. The distinct density dependences are attributed to the different impacts of the sp$^2$ and sp$^3$ motifs on the spectral heat capacity, vibrational mean free paths and group velocity. We additionally highlight the significant role of structural order in regulating the thermal conductivity of disordered carbon.

The glassy thermal conductivities observed in crystalline inorganic perovskites such as Cs$_3$Bi$_2$I$_6$Cl$_3$ is perplexing and lacking theoretical explanations. Here, we first experimentally measure such its thermal transport behavior from 20~K to 300~K, after synthesizing Cs$_3$Bi$_2$I$_6$Cl$_3$ single crystals. Using path-integral molecular dynamics simulations driven by machine learning potentials, we reveal that Cs$_3$Bi$_2$I$_6$Cl$_3$ has large lattice distortions at low temperatures, which may be related to the large atomic size mismatch. Employing the regularized Wigner thermal transport theory, we reproduce the experimental thermal conductivities based on lattice-distorted structures. This study thus provides a framework for predicting and understanding glassy thermal transport in materials with strong lattice disorder.

One of the most effective methods to enhance the strength of aluminum alloys involves modifying grain boundaries (GBs) through solute segregation. However, the fundamental mechanisms of solute segregation and their impacts on material properties remain elusive. In this study, we implemented highly efficient hybrid Monte Carlo and molecular dynamics (MCMD) algorithms in the graphics process units molecular dynamics (GPUMD) package. Using this efficient MCMD approach combined with a general-purpose machine-learning-based neuroevolution potential (NEP) for 16 elemental metals and their alloys, we simulated the segregation of 15 solutes in polycrystalline Al. Our results elucidate the segregation behavior and trends of 15 solutes in polycrystalline Al. Additionally, we investigated the impact of solutes on the strength of polycrystalline Al. The mechanisms underlying solute strengthening and embrittlement were analyzed at the atomistic level, revealing the importance of GB cohesion, as well as the nucleation and movement of Shockley dislocations, in determining the material's strength. We anticipate that our developed methods, along with our insights into solute segregation behavior in polycrystalline Al, will be valuable for the design of Al alloys and other multi-component materials, including medium-entropy materials, high-entropy materials, and complex concentrated alloys., Comment: 10 pages, 6 figures

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., Comment: 25 pages, 9 figures. This paper is part of the special topic, Machine Learning for Thermal Transport

Machine learned potentials (MLPs) have been widely employed in molecular dynamics (MD) simulations to study thermal transport. However, literature results indicate that MLPs generally underestimate the lattice thermal conductivity (LTC) of typical solids. Here, we quantitatively analyze this underestimation in the context of the neuroevolution potential (NEP), which is a representative MLP that balances efficiency and accuracy. Taking crystalline silicon, GaAs, graphene, and PbTe as examples, we reveal that the fitting errors in the machine-learned forces against the reference ones are responsible for the underestimated LTC as they constitute external perturbations to the interatomic forces. Since the force errors of a NEP model and the random forces in the Langevin thermostat both follow a Gaussian distribution, we propose an approach to correcting the LTC by intentionally introducing different levels of force noises via the Langevin thermostat and then extrapolating to the limit of zero force error. Excellent agreement with experiments is obtained by using this correction for all the prototypical materials over a wide range of temperatures. Based on spectral analyses, we find that the LTC underestimation mainly arises from increased phonon scatterings in the low-frequency region caused by the random force errors.

Infrared and Raman spectroscopy are widely used for the characterization of gases, liquids, and solids, as the spectra contain a wealth of information concerning in particular the dynamics of these systems. Atomic scale simulations can be used to predict such spectra but are often severely limited due to high computational cost or the need for strong approximations that limit application range and reliability. Here, we introduce a machine learning (ML) accelerated approach that addresses these shortcomings and provides a significant performance boost in terms of data and computational efficiency compared to earlier ML schemes. To this end, we generalize the neuroevolution potential approach to enable the prediction of rank one and two tensors to obtain the tensorial neuroevolution potential (TNEP) scheme. We apply the resulting framework to construct models for the dipole moment, polarizability, and susceptibility of molecules, liquids, and solids, and show that our approach compares favorably with several ML models from the literature with respect to accuracy and computational efficiency. Finally, we demonstrate the application of the TNEP approach to the prediction of infrared and Raman spectra of liquid water, a molecule (PTAF-), and a prototypical perovskite with strong anharmonicity (BaZrO3). The TNEP approach is implemented in the free and open source software package GPUMD, which makes this methodology readily available to the scientific community., Comment: 11 pages, 6 figures

Machine-learned potentials (MLPs) have exhibited remarkable accuracy, yet the lack of general-purpose MLPs for a broad spectrum of elements and their alloys limits their applicability. Here, we present a feasible approach for constructing a unified general-purpose MLP for numerous elements, demonstrated through a model (UNEP-v1) for 16 elemental metals and their alloys. To achieve a complete representation of the chemical space, we show, via principal component analysis and diverse test datasets, that employing one-component and two-component systems suffices. Our unified UNEP-v1 model exhibits superior performance across various physical properties compared to a widely used embedded-atom method potential, while maintaining remarkable efficiency. We demonstrate our approach's effectiveness through reproducing experimentally observed chemical order and stable phases, and large-scale simulations of plasticity and primary radiation damage in MoTaVW alloys. This work represents a significant leap towards a unified general-purpose MLP encompassing the periodic table, with profound implications for materials science., Comment: Main text with 17 pages and 8 figures; supplementary with 26 figures and 4 tables; source code and training/test data available

The lattice thermal conductivity (LTC) of Ga$_2$O$_3$ is an important property due to the challenge in the thermal management of high-power devices. We develop machine-learned neuroevolution potentials for single-crystalline $\beta$-Ga$_2$O$_3$ and $\kappa$-Ga$_2$O$_3$, and apply them to perform homogeneous nonequilibrium molecular dynamics simulations to predict their LTCs. The LTC of $\beta$-Ga$_2$O$_3$ was determined to be 10.3 $\pm$ 0.2 W/(m K), 19.9 $\pm$ 0.2 W/(m K), and 12.6 $\pm$ 0.2 W/(m K) along [100], [010], and [001], respectively, aligning with previous experimental measurements. For the first time, we predict the LTC of $\kappa$-Ga$_2$O$_3$ along [100], [010], and [001] to be 4.5 $\pm$ 0.0 W/(m K), 3.9 $\pm$ 0.0 W/(m K), and 4.0 $\pm$ 0.1 W/(m K), respectively, showing a nearly isotropic thermal transport property. The reduced LTC of $\kappa$-Ga$_2$O$_3$ versus $\beta$-Ga$_2$O$_3$ stems from its restricted low-frequency phonons up to 5 THz. Furthermore, we find that the $\beta$ phase exhibits a typical temperature dependence slightly stronger than $\sim T^{-1}$, whereas the $\kappa$ phase shows a weaker temperature dependence, ranging from $\sim T^{-0.5}$ to $\sim T^{-0.7}$., Comment: 8 pages, 7 figures

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 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 thermoelectric transport properties of a graphene antidot lattice., Comment: 8 pages, 4 figures

Amorphous silica (a-SiO$_2$) is a foundational disordered material for which the thermal transport properties are important for various applications. To accurately model the interatomic interactions in classical molecular dynamics (MD) simulations of thermal transport in a-SiO$_2$, we herein develop an accurate yet highly efficient machine-learned potential model that allowed us to generate a-SiO$_2$ samples closely resembling experimentally produced ones. Using the homogeneous nonequilibrium MD method and a proper quantum-statistical correction to the classical MD results, quantitative agreement with experiments is achieved for the thermal conductivities of bulk and 190 nm-thick a-SiO$_2$ films over a wide range of temperatures. To interrogate the thermal vibrations at different temperatures, we calculated the current correlation functions corresponding to the transverse acoustic (TA) and longitudinal acoustic (LA) collective vibrations. The results reveal that below the Ioffe-Regel crossover frequency, phonons as well-defined excitations, remain applicable in a-SiO$_2$ and play a predominant role at low temperatures, resulting in a temperature-dependent increase in thermal conductivity. In the high-temperature region, more phonons are excited, accompanied by a more intense liquid-like diffusion event. We attribute the temperature-independent thermal conductivity in the high-temperature range of a-SiO$_2$ to the collaborative involvement of excited phonon scattering and liquid-like diffusion in heat conduction. These findings provide physical insights into the thermal transport of a-SiO$_2$ and are expected to be applied to a vast range of amorphous materials., Comment: 12 pages, 7 figures in main text; 15 pages, 12 figures in Supplemental Material

Machine-learned potentials (MLPs) have become a popular approach of modelling interatomic interactions in atomistic simulations, but to keep the computational cost under control, a relatively short cutoff must be imposed, which put serious restrictions on the capability of the MLPs for modelling relatively long-ranged dispersion interactions. In this paper, we propose to combine the neuroevolution potential (NEP) with the popular D3 correction to achieve a unified NEP-D3 model that can simultaneously model relatively short-ranged bonded interactions and relatively long-ranged dispersion interactions. We show the improved descriptions of the binding and sliding energies in bilayer graphene can be obtained by the NEP-D3 approach compared to the pure NEP approach. We implement the D3 part into the GPUMD package such that it can be used out of the box for many exchange-correlation functionals. As a realistic application, we show that dispersion interactions result in approximately a 10% reduction in thermal conductivity for three typical metal-organic frameworks., Comment: 7 pages, 5 figures

Materials with low thermal conductivity usually have complex crystal structures. Herein we experimentally find that a simple crystal structure material AgTlI2 (I4/mcm) owns an extremely low thermal conductivity of 0.25 W/mK at room temperature. To understand this anomaly, we perform in-depth theoretical studies based on ab initio molecular dynamics simulations and anharmonic lattice dynamics. We find that the unique atomic arrangement and weak chemical bonding provide a permissive environment for strong oscillations of Ag atoms, leading to a considerable rattling behavior and giant lattice anharmonicity. This feature is also verified by the experimental probability density function refinement of single-crystal diffraction. The particularly strong anharmonicity breaks down the conventional phonon gas model, giving rise to non-negligible wavelike phonon behaviors in AgTlI2 at 300 K. Intriguingly, unlike many strongly anharmonic materials where a small propagative thermal conductivity is often accompanied by a large diffusive thermal conductivity, we find an unusual coexistence of ultralow propagative and diffusive thermal conductivities in AgTlI2 based on the thermal transport unified theory. This study underscores the potential of simple crystal structures in achieving low thermal conductivity and encourages further experimental research to enrich the family of materials with ultralow thermal conductivity., Comment: 6 pages 4 figures

The electron-phonon Wannier interpolation (EPWI) method is an efficient way to compute the properties of electron-phonon interactions (EPIs) accurately. This study presents a GPU-accelerated implementation of the EPWI method for computing transport properties, followed by a performance analysis. The implementation is tested on common systems such as aluminum and silicon. The results show complete consistency with those obtained through CPU computations. The proposed algorithm has the capability of computing the conductivity of aluminum in 20 minutes on a single NVIDIA Tesla V100 GPU, adopting a $200^3$ electron and phonon sampling grid. This speed is 173 times higher than the CPU-based algorithm, running on two nodes of the Intel Xeon Platinum 8260 CPU. Such impressive acceleration is achieved by carefully designing the algorithm to exploit the GPU's specific features. Furthermore, this methodology establishes a generic foundation for EPWI algorithms, which can be applied to other EPI-related properties.

Metal-organic frameworks (MOFs) are a family of materials that have high porosity and structural tunability and hold great potential in various applications, many of which requiring a proper understanding of the thermal transport properties. Molecular dynamics (MD) simulations play an important role in characterizing the thermal transport properties of various materials. However, due to the complexity of the structures, it is difficult to construct accurate empirical interatomic potentials for reliable MD simulations of MOFs. To this end, we develop a set of accurate yet highly efficient machine-learned potentials for three typical MOFs, including MOF-5, HKUST-1, and ZIF-8, using the neuroevolution potential approach as implemented in the GPUMD package, and perform extensive MD simulations to study thermal transport in the three MOFs. Although the lattice thermal conductivity (LTC) values of the three MOFs are all predicted to be smaller than 1 $\rm{W/(m\ K)}$ at room temperature, the phonon mean free paths (MFPs) are found to reach the sub-micrometer scale in the low-frequency region. As a consequence, the apparent LTC only converges to the diffusive limit for micrometer single crystals, which means that the LTC is heavily reduced in nanocrystalline MOFs. The sub-micrometer phonon MFPs are also found to be correlated with a moderate temperature dependence of LTC between those in typical crystalline and amorphous materials. Both the large phonon MFPs and the moderate temperature dependence of LTC fundamentally change our understanding of thermal transport in MOFs., Comment: 12 pages, 9 figures

Simulating collision cascades and radiation damage poses a long-standing challenge for existing interatomic potentials, both in terms of accuracy and efficiency. Machine-learning based interatomic potentials have shown sufficiently high accuracy for radiation damage simulations, but most existing ones are still not efficient enough to model high-energy collision cascades with sufficiently large space and time scales. To this end, we here extend the highly efficient neuroevolution potential (NEP) framework by combining it with the Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion potential, obtaining a NEP-ZBL framework. We train a NEP-ZBL model for tungsten and demonstrate its accuracy in terms of the elastic properties, melting point, and various energetics of defects that are relevant for radiation damage. We then perform large-scale molecular dynamics simulations with the NEP-ZBL model with up to 8.1 million atoms and 240 ps (using a single 40-GB A100 GPU) to study the difference of primary radiation damage in both bulk and thin-foil tungsten. While our findings for bulk tungsten are consistent with existing results simulated by embedded atom method (EAM) models, the radiation damage differs significantly in foils and shows that larger and more vacancy clusters as well as smaller and fewer interstitial clusters are produced due to the presence of a free surface.

It has recently been demonstrated that MoS2 with irregular interlayer rotations can achieve an extreme anisotropy in the lattice thermal conductivity (LTC), which is for example of interest for applications in waste heat management in integrated circuits. Here, we show by atomic scale simulations based on machine-learned potentials that this principle extends to other two-dimensional materials including C and BN. In all three materials introducing rotational disorder drives the through-plane LTC to the glass limit, while the in-plane LTC remains almost unchanged compared to the ideal bulk materials. We demonstrate that the ultralow through-plane LTC is connected to the collapse of their transverse acoustic modes in the through-plane direction. Furthermore, we find that the twist angle in periodic moir\'e structures representing rotational order provides an efficient means for tuning the through-plane LTC that operates for all chemistries considered here. The minimal through-plane LTC is obtained for angles between 1 and 4 degree depending on the material, with the biggest effect in MoS2. The angular dependence is correlated with the degree of stacking disorder in the materials, which in turn is connected to the slip surface. This provides a simple descriptor for predicting the optimal conditions at which the LTC is expected to become minimal., Comment: 10 pages, 6 figures

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., Comment: 8 pages, 7 figures

Polarizability is closely related to many fundamental characteristics of molecular systems and plays an indispensable role in simulating the Raman spectra. However, the calculations of polarizability for large systems still suffers from the limitations of processing ability of the quantum mechanical (QM) methods. This work assessed and compared the accuracy of the bond polarizability model (BPM) and a ML-based atomic polarizability model (AlphaML) in predicting polarizability of alkanes and then also investigated the ability of simulating Raman spectra. We found that the AlphaML has appreciable advantages over the BPM in learning the polarizability in the training data set and predicting polarizability of molecules that configurational differently from training structures. In addition, the BPM has inherent disadvantages in predicting polarizability anisotropy due to many factors including large uncertainties of estimating bond anisotropy, omitting of off-diagonal parameters in the construction of the model. As a result, the BPM has larger errors than the AlphaML in the simulation of anisotropic Raman scattering. Finally, we demonstrated that both the BPM and AlphaML suffer from transference to alkanes larger than those used in the training data sets, but the problem for the AlphaML can be circumvented by exploring more proper training structures., Comment: 37 pages, 10 figures

We report a local minimum in thermal conductivity in twisted bilayer graphene (TBG) at the angle of 1.08$^\circ$, which corresponds to the 'magic angle' in the transition of several other reported properties. Within the supercell of a moir\'e lattice, different stacking modes generate phonon scattering sites which reduce the thermal conductivity of TBG. The thermal magic angle arises from the competition between the delocalization of atomic vibrational amplitudes and stresses on one hand, and the increased AA stacking density on the other hand. The former effect weakens the scattering strength of a single scatterer while the latter one increases the density of scatterers. The combination of these two effects eventually leads to the apparition of the highlighted irregularity in heat conduction. The manifestation of a magic angle, disclosing new thermal mechanisms at nanoscale, further uncovers the unique physics of two-dimensional materials., Comment: 15 pages, 5 figures

Recently a novel two-dimensional (2D) C$_{60}$ based crystal called quasi-hexagonal-phase fullerene (QHPF) has been fabricated and demonstrated to be a promising candidate for 2D electronic devices [Hou et al. Nature 606, 507-510 (2022)]. We construct an accurate and transferable machine-learned potential to study heat transport and related properties of this material, with a comparison to the face-centered-cubic bulk-phase fullerene (BPF). Using the homogeneous nonequilibrium molecular dynamics and the related spectral decomposition methods, we show that the thermal conductivity in QHPF is anisotropic, which is 137(7) W/mK at 300 K in the direction parallel to the cycloaddition bonds and 102(3) W/mK in the perpendicular in-plane direction. By contrast, the thermal conductivity in BPF is isotropic and is only 0.45(5) W/mK. We show that the inter-molecular covalent bonding in QHPF plays a crucial role in enhancing the thermal conductivity in QHPF as compared to that in BPF. The heat transport properties as characterized in this work will be useful for the application of QHPF as novel 2D electronic devices., Comment: 11 pages, 12 figures

In a previous paper [Physical Review B \textbf{103}, 035417 (2021)], we showed that the equilibrium molecular dynamics (EMD) method can be used to compute the apparent thermal conductivity of finite systems. It has been shown that the apparent thermal conductivity from EMD for a system with domain length $2L$ is equal to that from nonequilibrium molecular dynamics (NEMD) for a system with domain length $L$. Taking monolayer silicence with an accurate machine learning potential as an example, here we show that the thermal conductivity values from EMD and NEMD agree for the same domain length if the NEMD is applied with periodic boundary conditions in the transport direction. Our results thus establish an exact equivalence between EMD and NEMD for thermal conductivity calculations., Comment: 5 pages, 5 figures

Phosphorus has diverse chemical bonds and even in its two-dimensional form there are three stable allotropes: black phosphorene (Black-P), blue phosphorene (Blue-P), and violet phosphorene (Violet-P). Due to the complexity of these structures, no efficient and accurate classical interatomic potential has been developed for them. In this paper, we develop an efficient machine-learned neuroevolution potential model for these allotropes and apply it to study thermal transport in them via extensive molecular dynamics (MD) simulations. Based on the homogeneous nonequilibrium MD method, the thermal conductivities are predicted to be $12.5 \pm 0.2$ (Black-P in armchair direction), $78.4 \pm 0.4$ (Black-P in zigzag direction), $128 \pm 3$ (Blue-P), and $2.36 \pm 0.05$ (Violet-P) $\mathrm{Wm^{-1}K^{-1}}$. The underlying reasons for the significantly different thermal conductivity values in these allotropes are unraveled through spectral decomposition, phonon eigenmodes, and phonon participation ratio. Under external tensile strain, the thermal conductivity in black-P and violet-P are finite, while that in blue-P appears unbounded due to the linearization of the flexural phonon dispersion that increases the phonon mean free paths in the zero-frequency limit., Comment: 10 pages, 10 figures, code and data available

Amorphous silicon (a-Si) is an important thermal-management material and also serves as an ideal playground for studying heat transport in strongly disordered materials. Theoretical prediction of the thermal conductivity of a-Si in a wide range of temperatures and sample sizes is still a challenge. Herein we present a systematic investigation of the thermal transport properties of a-Si by employing large-scale molecular dynamics (MD) simulations with an accurate and efficient machine-learned neuroevolution potential (NEP) trained against abundant reference data calculated at the quantum-mechanical density-functional-theory level. The high efficiency of NEP allows us to study the effects of finite size and quenching rate in the formation of a-Si in great detail. We find that it requires a simulation cell up to $64,000$ atoms (a cubic cell with a linear size of 11 nm) and a quenching rate down to $10^{11}$ K s$^{-1}$ for fully convergent thermal conductivity. Structural properties, including short- and medium-range order as characterized by the pair correlation function, angular distribution function, coordination number, ring statistics and structure factor are studied to demonstrate the accuracy of NEP and to further evaluate the role of quenching rate. Using both the heterogeneous and the homogeneous nonequilibrium MD methods and the related spectral decomposition techniques, we calculate the temperature- and thickness-dependent thermal conductivity values of a-Si and show that they agree well with available experimental results from 10 K to room temperature. Our results also highlight the importance of quantum effects in the calculated thermal conductivity and support the quantum correction method based on the spectral thermal conductivity.

We present our latest advancements of machine-learned potentials (MLPs) based on the neuroevolution potential (NEP) framework introduced in [Fan et al., Phys. Rev. B 104, 104309 (2021)] and their implementation in the open-source package GPUMD. We increase the accuracy of NEP models both by improving the radial functions in the atomic-environment descriptor using a linear combination of Chebyshev basis functions and by extending the angular descriptor with some four-body and five-body contributions as in the atomic cluster expansion approach. We also detail our efficient implementation of the NEP approach in graphics processing units as well as our workflow for the construction of NEP models, and we demonstrate their application in large-scale atomistic simulations. By comparing to state-of-the-art MLPs, we show that the NEP approach not only achieves above-average accuracy but also is far more computationally efficient. These results demonstrate that the GPUMD package is a promising tool for solving challenging problems requiring highly accurate, large-scale atomistic simulations. To enable the construction of MLPs using a minimal training set, we propose an active-learning scheme based on the latent space of a pre-trained NEP model. Finally, we introduce three separate Python packages, GPYUMD, CALORINE, and PYNEP, which enable the integration of GPUMD into Python workflows., Comment: 29 pages, 15 figures, code and data available

We study the interfacial thermal conductance of grain boundaries (GBs) between monolayer graphene and hexagonal boron nitride (h-BN) sheets using a combined atomistic approach. First, realistic samples containing graphene/h-BN GBs with different tilt angles are generated using the phase-field crystal (PFC) model developed recently [P. Hirvonen \textit{et al.}, Phys. Rev. B \textbf{100}, 165412 (2019)] that captures slow diffusive relaxation inaccessible to molecular dynamics (MD) simulations. Then, large-scale MD simulations using the efficient GPUMD package are performed to assess heat transport and rectification properties across the GBs. We find that lattice mismatch between the graphene and h-BN sheets plays a less important role in determining the interfacial thermal conductance as compared to the tilt angle. In addition, we find no significant thermal rectification effects for these GBs., Comment: 8 pages, 8 figures

In a previous paper [Fan Z \textit{et al}. 2021 Phys. Rev. B, \textbf{104}, 104309], we developed the neuroevolution potential (NEP), a framework of training neural network based machine-learning potentials using a natural evolution strategy and performing molecular dynamics (MD) simulations using the trained potentials. The atom-environment descriptor in NEP was constructed based on a set of radial and angular functions. For multi-component systems, all the radial functions between two atoms are multiplied by some fixed factors that depend on the types of the two atoms only. In this paper, we introduce an improved descriptor for multi-component systems, in which different radial functions are multiplied by different factors that are also optimized during the training process, and show that it can significantly improve the regression accuracy without increasing the computational cost in MD simulations., Comment: 7 pages, 8 figures, code and data available

We study the structural and mechanical properties of nanoporous (NP) carbon materials by extensive atomistic machine-learning (ML) driven molecular dynamics (MD) simulations. To this end, we retrain a ML Gaussian approximation potential (GAP) for carbon by recalculating the a-C structural database of Deringer and Cs\'anyi [Phys. Rev. B 2017, 95, 094203] adding van der Waals interactions. Our GAP enables a notable speedup and improves the accuracy of energy and force predictions. We use the GAP to thoroughly study the atomistic structure and pore-size distribution in computational NP carbon samples. These samples are generated by a melt-graphitization-quench MD procedure over a wide range of densities (from 0.5 to 1.7 g/cm$^3$) with structures containing 131,072 atoms. Our results are in good agreement with experimental data for the available observables, and provide a comprehensive account of structural (radial and angular distribution functions, motif and ring counts, X-ray diffraction patterns, pore characterization) and mechanical (elastic moduli and their evolution with density) properties. Our results show relatively narrow pore-size distributions, where the peak position and width of the distributions are dictated by the mass density of the materials. Our data allow further work on computational characterization of NP carbon materials, in particular for energy-storage applications, as well as suggest future experimental characterization of NP carbon-based materials.

We develop a neuroevolution-potential (NEP) framework for generating neural network based machine-learning potentials. They are trained using an evolutionary strategy for performing large-scale molecular dynamics (MD) simulations. A descriptor of the atomic environment is constructed based on Chebyshev and Legendre polynomials. The method is implemented in graphic processing units within the open-source GPUMD package, which can attain a computational speed over $10^7$ atom-step per second using one Nvidia Tesla V100. Furthermore, per-atom heat current is available in NEP, which paves the way for efficient and accurate MD simulations of heat transport in materials with strong phonon anharmonicity or spatial disorder, which usually cannot be accurately treated either with traditional empirical potentials or with perturbative methods., Comment: 18 pages, 10 figures, 2 tables, code and data available

Anomalous heat transport in one-dimensional nanostructures, such as nanotubes and nanowires, is a widely debated problem in condensed matter and statistical physics, with contradicting pieces of evidence from experiments and simulations. Using a comprehensive modeling approach, comprised of lattice dynamics and molecular dynamics simulations, we proved that the infinite length limit of the thermal conductivity of a (10,0) single-wall carbon nanotube is finite but this limit is reached only for macroscopic lengths due to thermal phonon mean free path of several millimeters. Our calculations showed that the extremely high thermal conductivity of this system at room temperature is dictated by quantum effects. Modal analysis showed that the divergent nature of thermal conductivity, observed in one-dimensional model systems, is suppressed in carbon nanotubes by anharmonic scattering channels provided by the flexural and optical modes with polarization in the plane orthogonal to the transport direction., Comment: 7 pages, 4 figures

The design of new applications, especially those based on heterogeneous integration, must rely on detailed knowledge of material properties, such as thermal conductivity (TC). To this end, multiple methods have been developed to study TC as a function of vibrational frequency. Here, we compare three spectral TC methods based on velocity decomposition in homogenous molecular dynamics simulations: Green-Kubo modal analysis (GKMA), the spectral heat current (SHC) method, and a method we propose called homogeneous nonequilibrium modal analysis (HNEMA). First, we derive a convenient per-atom virial expression for systems described by general many-body potentials, enabling compact representations of the heat current, each velocity decomposition method, and other related quantities. Next, we evaluate each method by calculating the spectral TC for carbon nanotubes, graphene, and silicon. We show that each method qualitatively agrees except at optical phonon frequencies, where a combination of mismatched eigenvectors and a large density of states produces artificial TC peaks for modal analysis methods. Our calculations also show that the HNEMA and SHC methods converge much faster than the GKMA method, with the SHC method being the most computationally efficient. Finally, we demonstrate that our single-GPU modal analysis implementation in GPUMD (Graphics Processing Units Molecular Dynamics) is over 1000 times faster than the existing LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) implementation on one CPU., Comment: 19 Pages, 8 Figures

We propose a way to properly interpret the apparent thermal conductivity obtained for finite systems using equilibrium molecular dynamics simulations (EMD) with fixed or open boundary conditions in the transport direction. In such systems the heat current autocorrelation function develops negative values after a correlation time which is proportional to the length of the simulation cell in the transport direction. Accordingly, the running thermal conductivity develops a maximum value at the same correlation time and eventually decays to zero. By comparing EMD with nonequilibrium molecular dynamics (NEMD) simulations, we conclude that the maximum thermal conductivity from EMD in a system with domain length 2L is equal to the thermal conductivity from NEMD in a system with domain length L. This facilitates the use of nonperiodic-boundary EMD for thermal transport in finite samples in close correspondence to NEMD., Comment: 7pages, 8 figures

Single-walled carbon nanotubes (SWCNTs) in their pristine form have high thermal conductivity whose further improvement has attracted a lot of interest. Some theoretical studies have suggested that the thermal conductivity of a $(10,10)$ SWCNT is dramatically enhanced by C$_{60}$ fullerene encapsulation. However, recent experiments on SWCNT bundles show that fullerene encapsulation leads to a reduction rather than an increase in thermal conductivity. Here, we employ three different molecular dynamics methods to study the influence of C$_{60}$ encapsulation on heat transport in a $(10,10)$ SWCNT. All the three methods consistently predict a reduction of the thermal conductivity of $(10,10)$ SWCNT upon C$_{60}$ encapsulation by $20\%-30\%$, in agreement with experimental results on bundles of SWCNTs. We demonstrate that there is a simulation artifact in the Green-Kubo method which gives anomalously large thermal conductivity from artificial convection. Our results show that the C$_{60}$ molecules conduct little heat compared to the outer SWCNT and reduce the phonon mean free paths of the SWCNT by inducing extra phonon scattering. We also find that the thermal conductivity of a $(10,10)$ SWCNT monotonically decreases with increasing filling ratio of C$_{60}$ molecules., Comment: 26 pages, 10 figures

Nanoconfinement induces many intriguing non-Fourier heat conduction phenomena that have been extensively studied in recent years, such as the nonlinear temperature profile inside the devices, the temperature jumps near the contacts, and the finite-size effects. The understanding of these phenomena, however, has been a matter of debate over the past two decades. In this work, we demonstrate a unified phonon interpretation of non-Fourier heat conduction which can help to understand these phenomena by a mode-to-mode correspondence between the non-equilibrium molecular dynamics (NEMD) simulations and the mode-resolved phonon Boltzmann transport equation (BTE). It is found that the nanoscale phonon transport characteristics including temperature profile, the heat flux value and the modal temperature depend on the applied thermal reservoirs on the two contacts. Our NEMD simulations demonstrate that Langevin thermostat behaves like an infinitely large thermal reservoir and provides thermally equilibrium mode-resolved phonon outlets, while biased reservoirs, e.g., Nose-Hoover chain thermostat and velocity rescaling method behave like non-equilibrium phonon outlets. Our interpretation clearly demonstrates that the non-Fourier heat transport phenomena are originated from a combination of non-diffusive phonon transport and phonon thermal nonequilibrium. This work provides a clear understanding of nanoscale heat transport and may guide the measurement and control of thermal transport in various applications., Comment: 23 pages, 12 figures

Silicon is an important material and many empirical interatomic potentials have been developed for atomistic simulations of it. Among them, the Tersoff potential and its variants are the most popular ones. However, all the existing Tersoff-like potentials fail to reproduce the experimentally measured thermal conductivity of diamond silicon. Here we propose a modified Tersoff potential and develop an efficient open source code called GPUGA (graphics processing units genetic algorithm) based on the genetic algorithm and use it to fit the potential parameters against energy, virial and force data from quantum density functional theory calculations. This potential, which is implemented in the efficient open source GPUMD (graphics processing units molecular dynamics) code, gives significantly improved descriptions of the thermal conductivity and phonon dispersion of diamond silicon as compared to previous Tersoff potentials and at the same time well reproduces the elastic constants. Furthermore, we find that quantum effects on the thermal conductivity of diamond silicon at room temperature are non-negligible but small: using classical statistics underestimates the thermal conductivity by about 10\% as compared to using quantum statistics., Comment: 9 pages, 6 figures

Atomically thin 2-dimensional heterostructures are a promising, novel class of materials with groundbreaking properties. The possiblity of choosing the many constituent components and their proportions allows optimizing these materials to specific requirements. The wide adaptability comes with a cost of large parameter space making it hard to experimentally test all the possibilities. Instead, efficient computational modelling is needed. However, large range of relevant time and length scales related to physics of polycrystalline materials poses a challenge for computational studies. To this end, we present an efficient and flexible phase-field crystal model to describe the atomic configurations of multiple atomic species and phases coexisting in the same physical domain. We extensively benchmark the model for two-dimensional binary systems in terms of their elastic properties and phase boundary configurations and their energetics. As a concrete example, we demonstrate modelling lateral heterostructures of graphene and hexagonal boron nitride. We consider both idealized bicrystals and large-scale systems with random phase distributions. We find consistent relative elastic moduli and lattice constants, as well as realistic continuous interfaces and faceted crystal shapes. Zigzag-oriented interfaces are observed to display the lowest formation energy., Comment: 17 pages, 13 figures

Nonequilibrium molecular dynamics (NEMD) has been extensively used to study thermal transport at various length scales in many materials. In this method, two local thermostats at different temperatures are used to generate a nonequilibrium steady state with a constant heat flux. Conventionally, the thermal conductivity of a finite system is calculated as the ratio between the heat flux and the temperature gradient extracted from the linear part of the temperature profile away from the local thermostats. Here we show that, with a proper choice of the thermostat, the nonlinear part of the temperature profile should actually not be excluded in thermal transport calculations. We compare NEMD results against those from the atomistic Green's function method in the ballistic regime, and those from the homogeneous nonequilibrium molecular dynamics method in the ballistic-to-diffusive regime. These comparisons suggest that in all the transport regimes, one should directly calculate the thermal conductance from the temperature difference between the heat source and sink and, if needed, convert it to the thermal conductivity by multiplying it with the system length. Furthermore, we find that the Langevin thermostat outperforms the Nos\'{e}-Hoover (chain) thermostat in NEMD simulations because of its stochastic and local nature. We show that this is particularly important for studying asymmetric carbon-based nanostructures, for which the Nos\'{e}-Hoover thermostat can produce artifacts leading to unphysical thermal rectification. Our findings are important to obtain correct results from molecular dynamics simulations of nanoscale heat transport as the accuracy of the interatomic potentials is rapidly improving., Comment: 13 pages, 12 figures

Understanding the mechanisms of thermal conduction in graphene is a long-lasting research topic, due to its high thermal conductivity. Peierls-Boltzmann transport equation (PBTE) based studies have revealed many unique phonon transport properties in graphene, but most previous works only considered three-phonon scatterings and relied on interatomic force constants (IFCs) extracted at 0 K. In this paper, we explore the roles of four-phonon scatterings and the temperature dependent IFCs on phonon transport in graphene through our PBTE calculations. We demonstrate that the strength of four-phonon scatterings would be severely overestimated by using the IFCs extracted at 0 K compared with those corresponding to a finite temperature, and four-phonon scatterings are found to significantly reduce the thermal conductivity of graphene even at room temperature. In order to reproduce the prediction from molecular dynamics simulations, phonon frequency broadening has to be taken into account when determining the phonon scattering rates. Our study elucidates the phonon transport properties of graphene at finite temperatures, and could be extended to other crystalline materials., Comment: 31 pages, 8 figures

In recent years, predictive computational modeling has become a cornerstone for the study of fundamental electronic, optical, and thermal properties in complex forms of condensed matter, including Dirac and topological materials. The simulation of quantum transport in realistic materials calls for the development of linear scaling, or order-$N$, numerical methods, which then become enabling tools for guiding experimental research and for supporting the interpretation of measurements. In this review, we describe and compare different order-$N$ computational methods that have been developed during the past twenty years, and which have been used extensively to explore quantum transport phenomena in disordered media. We place particular focus on the zero-frequency electrical conductivities derived within the Kubo-Greenwood and Kubo-Streda formalisms, and illustrate the capabilities of these methods to tackle the quasi-ballistic, diffusive, and localization regimes of quantum transport in the noninteracting limit. The fundamental issue of computational cost versus accuracy of various proposed numerical schemes is addressed in depth. We then illustrate the usefulness of these methods with various examples of transport in disordered materials, such as polycrystalline and defected graphene models, 3D metals and Dirac semimetals, carbon nanotubes, and organic semiconductors. Finally, we extend the review to the study of spin dynamics and topological transport, for which efficient approaches for calculating charge, spin, and valley Hall conductivities are described., Comment: 54 pages, 31 figures, Invited Review of Physics Reports

Thermal properties of molybdenum disulfide (MoS$_2$) have recently attracted attention related to fundamentals of heat propagation in strongly anisotropic materials, and in the context of potential applications to optoelectronics and thermoelectrics. Multiple empirical potentials have been developed for classical molecular dynamics (MD) simulations of this material, but it has been unclear which provides the most realistic results. Here, we calculate lattice thermal conductivity of single- and multi-layer pristine MoS$_2$ by employing three different thermal transport MD methods: equilibrium, nonequilibrium, and homogeneous nonequilibrium ones. These methods allow us to verify the consistency of our results and also facilitate comparisons with previous works, where different schemes have been adopted. Our results using variants of the Stillinger-Weber potential are at odds with some previous ones and we analyze the possible origins of the discrepancies in detail. We show that, among the potentials considered here, the reactive empirical bond order (REBO) potential gives the most reasonable predictions of thermal transport properties as compared to experimental data. With the REBO potential, we further find that isotope scattering has only a small effect on thermal conduction in MoS$_2$ and the in-plane thermal conductivity decreases with increasing layer number and saturates beyond about three layers. We identify the REBO potential as a transferable empirical potential for MD simulations of MoS$_2$ which can be used to study thermal transport properties in more complicated situations such as in systems containing defects or engineered nanoscale features. This work establishes a firm foundation for understanding heat transport properties of MoS$_2$ using MD simulations., Comment: 14 pages, 6 figures

We use a phase field crystal model to generate large-scale bicrystalline and polycrystalline single-layer hexagonal boron nitride (h-BN) samples and employ molecular dynamics (MD) simulations with the Tersoff many-body potential to study their heat transport properties. The Kapitza thermal resistance across individual h-BN grain boundaries is calculated using the inhomogeneous nonequilibrium MD method. The resistance displays strong dependence on the tilt angle, the line tension and the defect density of the grain boundaries. We also calculate the thermal conductivity of pristine h-BN and polycrystalline h-BN with different grain sizes using an efficient homogeneous nonequilibrium MD method. The in-plane and the out-of-plane (flexural) phonons exhibit different grain size scalings of the thermal conductivity in polycrystalline h-BN and the extracted Kapitza conductance is close to that of large-tilt-angle grain boundaries in bicrystals., Comment: 12 pages, 9 figures, 2 tables, bicrystalline and polycrystalline samples available

While the microscopic structure of defected solid crystalline materials has significant impact on their physical properties, efficient and accurate determination of a given polycrystalline microstructure remains a challenge. In this paper we present a highly generalizable and reliable variational method to achieve this goal for two-dimensional crystalline and quasicrystalline materials. The method is benchmarked and optimized successfully using a variety of large-scale systems of defected solids, including periodic structures and quasicrystalline symmetries to quantify their microstructural characteristics, e.g., grain size and lattice misorientation distributions. We find that many microstructural properties show universal features independent of the underlying symmetries.