Homogeneous Nonequilibrium Molecular Dynamics Research Articles

Thermal conductivity calculation using homogeneous non-equilibrium molecular dynamics simulation with Allegro

In this study, we derive the heat flux formula for the Allegro model, one of machine-learning interatomic potentials using the equivariant deep neural network, to calculate lattice thermal conductivity using the homogeneous non-equilibrium molecular dynamics (HNEMD) method based on the Green-Kubo formula. Allegro can construct more advanced atomic descriptors than conventional ones, and can be applied to multicomponent and large-scale systems, providing a significant advantage in estimating the thermal conductivity of anharmonic materials, such as thermoelectric materials. In addition, the spectral heat current (SHC) method, recently developed for the HNEMD framework (HNEMD-SHC), allows the calculation of not only the total thermal conductivity but also its frequency components. The verification of the heat flux and the demonstration of HNEMD-SHC method are performed for the extremely anharmonic low-temperature phase of Ag2Se.

Thermal transport in disordered wurtzite ScAlN alloys using machine learning interatomic potentials

Disordered wurtzite ScAlN alloy is an emerging material with a strong piezoelectric coefficient that promotes its application in 5 G communication devices. To extend the service lifetime of the device, it is an urgent task to understand the thermal transport properties of this alloy for thermal management. Yet a thorough understanding of the thermal transport in the disordered ScAlN alloy is hindered by the inaccuracy of empirical interatomic potentials adopted for atomistic simulations. Herein, we developed an accurate and efficient neuro-evolution potential and combined it with Allen and Feldman theory and molecular dynamics to investigate thermal transport properties of disordered ScAlN alloys with varying ratios. This potential is shown to have the same computational accuracy as first principles by comparing the phonon spectrum and radial distribution function. Based on the homogeneous non-equilibrium molecular dynamics simulations, the lattice thermal conductivity of the disordered alloys reaches its maximum value (less than 12 W/m·K) at room temperature. Further investigation of the temperature effects on the vibrational modes shows that it affects the mode lifetimes for the propagons and the mode heat capacities for the diffusions. In addition, it was found that both mode lifetime and diffusivity decrease with increasing Sc concentration. This work aims to train a universal machine learning interatomic potential of disordered ScAlN alloys, seeking to accelerate the thermal management design of related devices.

Thermal transports of 2D phosphorous carbides by machine learning molecular dynamics simulations

Carbon phosphide is a newly discovered two-dimensional semiconductor material which wrinkles and has a significant carrier mobility. Due to lack an accurate force field, the use of molecular dynamics to study its phonon-dominated thermal conductivity which lead to inaccurate results. At present, the use of machine learning to construct a high-precision force field has become the mainstream research method to solve this problem. The main work of this study is to construct a comprehensive training sets for Phosphorus-Doped Graphene (PCn) (n = 3, 5, 6) and to use the fitted potential to calculate the related thermal properties. The research found that (PC5) exhibited anisotropic behavior, with a thermal conductivity of 106.6 Wm−1K−1 in the y-direction and 63.6 Wm−1K−1 in the x-direction. In comparison, (PC6) and (PC3) showed isotropic behavior, with thermal conductivity of approximately 104 Wm−1K−1 and 76.83 Wm−1K−1, respectively. Compared to monolayer graphene, the lower thermal conductivity of PCn is mainly attributed to phonon-phonon scattering effects, which are limited by the regular wrinkled structure. Additionally, low-frequency phonon have been found to have a significant impact on the thermal performance of PCn. Furthermore, we investigated the influence of uniaxial strain on the PC6 and observed an increase in the thermal conductivity with increasing strain. This study used key computational and analytical techniques, including phonon dispersion relations, homogeneous nonequilibrium molecular dynamics method, spectral thermal conductivity analysis. These findings provide a theoretical basis for understanding the thermal transport properties of PCn and will guide its potential applications value.

Dissimilar thermal transport properties in κ-Ga2O3 and β-Ga2O3 revealed by homogeneous nonequilibrium molecular dynamics simulations using machine-learned potentials

The lattice thermal conductivity (LTC) of Ga2O3 is an important property due to the challenge in the thermal management of high-power devices. In this work, we develop machine-learned neuroevolution potentials (NEPs) for single-crystalline β−Ga2O3 and κ−Ga2O3 and demonstrate their accuracy in modeling thermal transport properties. Combining NEP-driven homogeneous non-equilibrium molecular dynamics simulations with tensor analysis, we determine the spatial distributions of LTCs for two Ga2O3 crystals, showing dissimilar thermal behaviors. Specifically, β−Ga2O3 shows isotropic thermal transport properties, with the LTCs along [100], [010], and [001] directions being predicted to be 10.3±0.2, 19.9±0.2, and 12.6±0.2 W/(m K), respectively, consistent with previous experimental measurements. For κ−Ga2O3, our predictions suggest nearly isotropic thermal transport properties, with the LTCs along [100], [010], and [001] being estimated to be 4.5±0.1, 3.9±0.1, and 4.0±0.1 W/(m K). The reduced LTC of κ−Ga2O3 vs β−Ga2O3 stems from its restricted low-frequency phonons up to 5 THz. Furthermore, we find that the β phase exhibits a typical temperature dependence slightly stronger than ∼T−1, whereas the κ phase shows a weaker temperature dependence, ranging from ∼T−0.5 to ∼T−0.7.

Accurate prediction of heat conductivity of water by a neuroevolution potential.

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.

Low lattice thermal conductivity with two-channel thermal transport in the superatomic crystal PH4AlBr4

Designing novel crystalline materials composed of light and nontoxic elements, but with ultralow thermal conductivity insensitive to temperature, has been a long-standing challenge. One effective strategy is to utilize superatoms as building blocks to introduce hierarchical bonding and vibration for suppressing phonon velocity and enhancing high-order scattering. However, far fewer theoretical efforts have been made to understand the lattice dynamics and thermal transport in superatom-based materials. Herein, different from most of the existing works reported hitherto on atom-superatom hybrid systems, we carry out a comprehensive study on the phonon interaction and thermal transport in superatomic crystal $\mathrm{P}{\mathrm{H}}_{4}\mathrm{Al}{\mathrm{Br}}_{4}$ solely consisting of superatoms (superalkali and superhalogen), by employing homogeneous nonequilibrium molecular dynamics with active-learning potential, combined with density functional theory and unified theory of phonon thermal transport. We find that the supersalt $\mathrm{P}{\mathrm{H}}_{4}\mathrm{Al}{\mathrm{Br}}_{4}$ crystal exhibits amorphouslike ultralow lattice thermal conductivity of $0.329\ensuremath{-}0.286\phantom{\rule{0.16em}{0ex}}\mathrm{W}\phantom{\rule{0.16em}{0ex}}{\mathrm{m}}^{\ensuremath{-}1}\phantom{\rule{0.16em}{0ex}}{\mathrm{K}}^{\ensuremath{-}1}$ from 200 to 600 K. In particular, due to the strong quartic anharmonicity, the contribution of four-phonon scattering reaches 38% of total scattering rates at 300 K, while the phonon coherence's contribution is comparable with that of population, resulting from significant phonon localization. These theoretical findings demonstrate that superatom-assembled materials exhibit features distinguished from atom-based crystals and provide a unique platform for exploring ultralow thermal conductivity with enhanced two-channel thermal transport.

Thermal conductivity across transition metal dichalcogenide bilayers

Cross-plane thermal conductivities of transition metal dichalcogenide (M= Mo, W; X= S, Se) bilayers are determined by homogeneous nonequilibrium molecular dynamics. The best insulator is found to be : , closely followed by : and : . Thermal conductivities of heterobilayers are close to the average of the two corresponding homobilayer values, despite the mass and lattice mismatch. To disentangle the effects of atomic mass, lattice constant, and interaction potential, these three parameters are systematically varied. Phonon spectral analysis provides further insight into their roles and reveals the weak influence of spectral overlap between the two layers and the dominance of boundary scattering. The observed trends can be rationalized using Slack's formula in terms of the average atomic mass and the Debye temperature. Accurate interlayer interaction potentials are developed based on experimental elastic constants. Their effect on the bilayer cross-plane thermal conductivities is found to be minor.

Anisotropic and high thermal conductivity in monolayer quasi-hexagonal fullerene: A comparative study against bulk phase fullerene

Recently a novel two-dimensional (2D) C60 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.

Thermal conductivity of fivefold twinned silicon-germanium heteronanowires.

The thermal transport properties of five-fold twinned (5FT) germanium-silicon (Ge-Si) heteronanowires (h-NWs) with varying cross-sectional areas, germanium (Ge) domain ratios and heterostructural patterns are investigated using homogeneous nonequilibrium molecular dynamics (HNEMD) simulations. The results demonstrate a distinctive behavior in the thermal conductivity (κ) of 5FT-NWs, characterized by a "flipped" trend at a critical cross-sectional area. This behavior is attributed to the hydrodynamic phonon flow, arising from the normal three-phonon scattering process in the low-frequency region. In addition, the composition ratio of 5FT-NWs has a significant impact on reducing the κ of 5FT-NWs and suppressing the hydrodynamic effect. Intriguingly, as the homogeneous element domains are separated, stronger phonon hydrodynamic flows are observed in comparison to the adjacent homogeneous element domains. By analyzing various phonon properties, including phonon dispersion, three-phonon scattering rate, and phonon mean free path, critical insights into the origin of the differential κ in different 5FT-NW structures are provided. The findings deepen the understanding of the thermal transport properties of nanomaterials and hold implications for the design and development of nanoelectronics and thermoelectric devices.

Isotope doping-induced crossover shift in the thermal conductivity of thin silicon nanowires

Here, using homogeneous nonequilibrium molecular dynamics simulations, we report the thermal transport characteristics of thin Si nanowires (NWs) with varying size and isotope doping ratio. It is identified that crossover in the thermal conductivity (κ) of both isotope doping-free and isotope doped Si-NWs appears at critical sizes, below which κ is enlarged with decreasing size because the hydrodynamic phonon flow predominates, above which, due to the dominant phonon boundary scattering, opposite behavior is observed. With increasing isotope doping, however, the critical size in minimizing the κ is moved to small values because the phonon impurity scattering caused by isotope doping is critically involved. Moreover, there is a critical isotope doping (<50%) in the critical size motion, originating from that, above which, the critical size no longer moves due to the persistence of hydrodynamic phonon flow. This study provides new insights into the thermal transport behaviors of quasi-1D structures.

Phonon thermal transport in graphene/h-BN superlattice monolayers

The superlattice monolayers composed of periodically assembled graphene and hexagonal boron nitride (h-BN) superlattices by strong covalent bonds have aroused great interest in academia and industry due to their adjustable physical properties, which are expected to play an important role in the modulation of material properties, especially in thermal transport properties. In general, thermal transport properties highly rely on the characteristic of nanostructures. Thus, understanding the impact of structural characteristics on physical properties is essential for the design of high-performance materials. In this study, we report our results from three different methods including non-equilibrium molecular dynamics, homogeneous non-equilibrium molecular dynamics, and spectral heat current decomposition to explore the phonon thermal transport properties in different graphene/h-BN superlattice monolayers, which are constructed by periodically stitched and equal-sized graphene and h-BN superlattices. We find that the thermal conductivities decrease first and then increase with the increase of periodic length, corresponding to a transition from coherent transport to incoherent transport, and a minimum thermal conductivity at a certain period is identified. These results provide a theoretical validation of a possible control on thermal properties at the nanoscale, which may have potential applications in thermoelectric devices once upon experimental validation.

Thermal transport of bilayer graphene: a homogeneous nonequilibrium molecular dynamics study

In this paper, the thermal transport of bilayer graphene is revisited by the homogeneous nonequilibrium molecular dynamics (HNEMD) method realized in a graphics processing unit based molecular dynamics package, GPUMD. Our simulations are carried out in three-dimensional boxes. An optimized Teroff potential is used to describe C-C covalent bonds, and the inter-layer van der Waals interaction is described by the 12-6 Lennard-Jones potential. Since the HNEMD method is homogeneous without boundary scattering, we obtain more accurate results than previous nonequilibrium molecular dynamics studies. Through intensive simulations, the main findings of the paper are as follows. (1) We find an efficient simulation setting, yielding results which are in good agreement with the experimental data; (2) Although the increase of the system size and the extension of production period has a beneficial effect to obtain more convergent results, these lead to large discrepancy with the experimental data; (3) The thermal conductivity is reduced by the phonon scattering among graphene layers; (4) The spectral decomposition of thermal conductivity shows that the bilayer graphene retains the spectrum of its monolayer counterpart except for the increasing magnitudes; (5) The spectral phonon mean free path of bilayer graphene reveals the presence of a new phonon mode which may have a negative contribution to the thermal conductivity.

Estimating thermal conductivity of α-Ag2Se using ANN potential with Chebyshev descriptor

We clarified computational conditions required to construct artificial neural network (ANN) potentials with a Chebyshev descriptor that accurately estimate thermal conductivity (TC). To this end, training data were generated by an empirical potential; α-Ag2Se was selected as the test system and TC was estimated by the homogeneous nonequilibrium molecular dynamics method based on the Green–Kubo formula. The ANN potentials trained with first-principles molecular dynamics data were constructed and the TCs were estimated based on the clarified conditions. Since ANN potentials for multicomponent systems are constructed efficiently with the Chebyshev descriptor, accurate TCs for such systems can be relatively easily estimated.

Equilibrium and nonequilibrium molecular dynamics methods to compute the first normal stress coefficient of a model polymer solution

The first normal pressure (or stress) difference is directly computed from the local values of the pressure tensor components in molecular dynamics simulations of planar Poiseuille flow for a low molecular weight polymeric fluid. The resulting zero shear rate normal pressure difference agrees very well with the value computed using homogeneous shear simulations and the SLLOD algorithm, and less well with the result of the Coleman-Markowitz equation evaluated at equilibrium. This resolves doubts about the effects of homogeneous thermostats in homogeneous nonequilibrium molecular dynamics simulations.

Thermal conductivity of monolayer hexagonal boron nitride: From defective to amorphous

Through extensive molecular dynamics simulations, we completed the thermal transport properties study of the monolayer hexagonal boron nitride (h-BN) films from the defect state to the amorphous state. To this end, a defective h-BN model construction program has been developed, which can customize the nature of the defects and realize the transition from the defect state to the amorphous state. By performing homogeneous non-equilibrium molecular dynamics simulations based on Tersoff multi-body potential, the thermal conductivity results of large-size h-BN films have been achieved. For the defect state h-BN, we studied the effect of different defect types and concentrations on its thermal conductivity and verified it from the phonon mechanism. After the transition to two-dimensional amorphous BN, we discussed and analyzed the influence of amorphous concentration, amorphous defect ratio and temperature on its thermal transport properties. On this basis, the restriction relationship between phonon-phonon scattering term and phonon-defect scattering term was also proposed from the atomic level mechanisms. Our results constitute a step in the deterministic engineering of thermal management devices in 2D materials, and hold great promise for the application of h-BN defect engineering in the field of nanomaterial thermal design.

Studies of weak external force on thermal conductivity of complex liquids using homogenous perturbed simulations

A modified Evan-Gillan homogenous perturbed scheme was employed to measure the thermal response and corresponding normalized plasma thermal conductivity (λ0) of strongly coupled dusty plasmas (SCDPs). Homogenous perturbed molecular dynamics (HPMD) approach was employed to analyze the efficiency and compared the obtained outcomes of perturbed heat flux vector to the outcomes computed through equilibrium molecular dynamics (EMD) approach. Our new outcomes show that the thermal response of heat energy current through HPMD and EMD approaches is an excellent agreement with each other for SCDPs at much low value of normalized external perturbation (F*). The obtained HPMD simulations demonstrate that the presented approach provides good outcomes with fast convergence and small system size for low-tointermediate plasma couplings (Γ) at very weak F*. It was found that the measured plasma λ0 outcomes at nearly equilibrium weak F* (= 0.00005) are in reasonable agreement with earlier results obtained from EMD, homogenous and inhomogeneous nonequilibrium molecular dynamics simulations and theoretical predictions and generally overpredicted the plasma λ0 by ~1%–~24% depending on the arrangement of plasma parameters (Г, к). It was shown that the HPMD is an excellent approach to calculate the plasma λ0 and to recognize the elementary patterns in 3D SCDPs.

Influence of thermostatting on nonequilibrium molecular dynamics simulations of heat conduction in solids.

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 into the thermal conductivity by multiplying it with the system length. Furthermore, we find that the Langevin thermostat outperforms the Nosé-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é-Hoover thermostat can produce artifacts leading to unphysical thermal rectification.

Crystalline order and disorder in dusty plasmas investigated by nonequilibrium molecular dynamics simulations**Project supported by the National Natural Science Foundation of China (Grant Nos. 11505286 and 11805272).

The particle structure of a complex system has been explored through a unique Evans's homogenous nonequilibrium molecular dynamics (HNEMD) simulation technique. The crystalline order–disorder structures (OD-structures) and the corresponding energies of three-dimensional (3D) nonideal complex systems (NICSs) have been measured over a wide range of plasma states (Γ, κ) for a body-centered cubic (BCC) structure. The projected technique provides accurate OD-structures with fast convergence and applicable to very small size effect for different temperatures ( 1 / and constant force field ( values. The OD-structure obtained through HNEMD approach is found to be reasonable agreement and more reliable than those earlier identified by simulation approaches and experimental data of NICSs. New simulations of OD-structures show that dusty plasma remains in crystalline (strongly coupled) state at lower temperature and constant values, for the whole simulation runs. Our investigations show that the crystalline structure is changed and the particle structure switches from intermediate to disorder (nonideal gaseous) state with an increase of the system's temperature. It has been shown that the long range order shifts toward lower temperature with increasing *. The presented technique exhibits that the potential energy has a maximum value when the dusty plasma remains in crystalline states (low temperatures), which confirms earlier 3D simulation results.

Homogeneous nonequilibrium molecular dynamics method for heat transport and spectral decomposition with many-body potentials

The standard equilibrium Green-Kubo and nonequilibrium molecular dynamics (MD) methods for computing thermal transport coefficients in solids typically require relatively long simulation times and large system sizes. To this end, we revisit here the homogeneous nonequilibrium MD method by Evans [Phys. Lett. A \textbf{91}, 457 (1982)] and generalize it to many-body potentials that are required for more realistic materials modeling. We also propose a method for obtaining spectral conductivity and phonon mean free path from the simulation data. This spectral decomposition method does not require lattice dynamics calculations and can find important applications in spatially complex structures. We benchmark the method by calculating thermal conductivities of three-dimensional silicon, two-dimensional graphene, and a quasi-one-dimensional carbon nanotube and show that the method is about one to two orders of magnitude more efficient than the Green-Kubo method. We apply the spectral decomposition method to examine the long-standing dispute over thermal conductivity convergence vs divergence in carbon nanotubes.

Thermal transport properties of single-layer black phosphorus from extensive molecular dynamics simulations

We compute the anisotropic in-plane thermal conductivity of suspended single-layer black phosphorus (SLBP) using three molecular dynamics (MD) based methods, including the equilibrium MD method, the nonequilibrium MD (NEMD) method, and the homogeneous NEMD (HNEMD) method. Two existing parameterizations of the Stillinger–Weber (SW) potential for SLBP are used. Consistent results are obtained for all the three methods and conflicting results from previous MD simulations are critically assessed. Among the three methods, the HNEMD method is the most and the NEMD method the least efficient. The thermal conductivity values from our MD simulations are about an order of magnitude larger than the most recent predictions obtained using the Boltzmann transport equation approach considering long-range interactions in density functional theory calculations, suggesting that the short-range SW potential might be inadequate for describing the phonon anharmonicity in SLBP.