Lennard-Jones Argon Research Articles

A constitutive relation generalizing the Navier–Stokes theory to high-rate regimes

We propose a constitutive equation for flows of gases in high-rate regimes where the Navier–Stokes theory breaks down. The model generalizes the Navier–Stokes relation and agrees well with that model in all lower rate flows examined. Our proposed constitutive relation is calibrated with the method of objective molecular dynamics (OMD) using families of compressible and incompressible flows of Lennard-Jones argon. The constitutive relation makes use of the higher-order objective strain rates due to Rivlin and Ericksen (J. Rat. Mech. Anal., vol. 4, 1955, pp. 323–425). The constitutive relation is fully frame-indifferent, and the macroscopic flows corresponding to the OMD simulations are exact solutions for the proposed model. The model is shown to agree with atomistic results much better than the Navier–Stokes equations in the transition regime. The success of our model indicates that it is not higher gradients that become important in the high-rate regime, but rather higher rates of change of the strain rate tensor. While somewhat more complicated to implement than the Navier–Stokes relation, the proposed model is expected to be compatible with existing methods of computational fluid dynamics and may extend those methods to higher rate regimes, while preserving their ability to handle large spatial scales.

Atomic dynamics in fluids: Normal mode analysis revisited.

Developing microscopic understanding of the thermal properties of liquids is challenging due to their strong dynamic disorder, which prevents characterization of the atomic degrees of freedom. There have been significant research interests in the past few decades to extend the normal mode analysis for solids to instantaneous structures of liquids. However, the nature of normal modes that arise from these unstable structures is still elusive. In this paper, we explore the instantaneous eigenmodes of dynamical matrices of various Lennard-Jones argon liquid and gas systems at high temperatures and show that the normal modes can be interpreted as an interpolation of T→∞ (gas) and T=0 (solid) mode descriptions. We find that normal modes become increasingly collisional and translational, recovering atomistic gaslike behavior rather than vibrational with increase in temperature, suggesting that normal modes in liquids may be described by both solidlike and gaslike modes.

Computing free energy barriers for the nucleation of complex network mesophases.

A previously introduced framework to identify local order parameters (OPs) distinctive of incipient complex mesophases, such as bicontinuous network phases, is used in this work to evaluate nucleation free-energy barriers. The sampling techniques considered are the mean-first-passage-time (MFPT) method and novel variants of umbrella sampling, including Hybrid Monte Carlo (HMC) and a dual-OP-method that uses a blunter global OP for the umbrella bias while keeping record of configurations for analysis with a local OP. These methods were chosen for their ability to minimize or avoid frequent calculation of the expensive local OP, which makes their continuous on-the-fly tracking computationally very inefficient. These techniques were first validated by studying phase-transition barriers of model systems, i.e., the vapor-liquid nucleation of Lennard-Jones argon and a binary nanoparticle model. The disorder-to-order free energy barrier was then traced for the double gyroid and single diamond formed by mesoscopic bead-spring macromolecular models. The dual OP method was found to be the most robust and computationally efficient, since, unlike HMC, it does not require the expensive local OP to be computed on-the-fly, and unlike the MFPT method, it can negotiate large barriers aided by the biased sampling. The dual OP method requires, however, that a cheap global OP be identified and correlated (in a post-processing step) with the local OP that describes the structure of the critical nucleus, a process that can be aided by machine learning.

Can we define a unique microscopic pressure in inhomogeneous fluids?

The estimation of a microscopic pressure tensor in an adsorbed thin film on a planar surface remains a challenge in both experiment and theory. While the normal pressure is well-defined for a planar surface, the tangential pressure at a point is not uniquely defined at the nanoscale. We report a new method that allows us to calculate the local pressure tensor and its spatial integral using an arbitrary contour definition of the "virial-route" local pressure tensor. We show that by integrating the local tangential pressure over a small region of space, roughly the range of the intermolecular forces, it is possible to define a coarse-grained tangential pressure that appears to be unique and free from ambiguities in the definition of the local pressure tensor. We support our argument by presenting the results for more than ten types of contour definitions of the local pressure tensor. By defining the coarse-grained tangential pressure, we can also find the effective thickness of the adsorbed layer and, in the case of a porous material, the statistical pore width. The coarse-grained in-layer and in-pore tangential pressures are determined for Lennard-Jones argon adsorbed in realistic carbon slit pores, providing a better understanding of the pressure enhancement for strongly wetting systems.

Phonon confinement and transport in ultrathin films

Thermal transport by phonons in films with thicknesses of less than 10 nm is investigated in a soft system (Lennard-Jones argon) and a stiff system (Tersoff silicon) using two-dimensional lattice dynamics calculations and the Boltzmann transport equation. This approach uses a unit cell that spans the film thickness, which removes approximations related to the finite cross-plane dimension required in typical three-dimensional-based approaches. Molecular dynamics simulations are performed to obtain finite-temperature structures for the lattice dynamics calculations and to predict thermal conductivity benchmarks. Thermal conductivity decreases with decreasing film thickness until the thickness reaches four unit cells (2.1 nm) for argon and three unit cells (1.6 nm) for silicon. With a further decrease in film thickness, thermal conductivity plateaus in argon while it increases in silicon. This unexpected behavior, which we identify as a signature of phonon confinement, is a result of an increased contribution from low-frequency phonons, whose density of states increases as the film thickness decreases. Phonon mode-level analysis suggests that confinement effects emerge below thicknesses of ten unit cells (5.3 nm) for argon and six unit cells (3.2 nm) for silicon. Thermal conductivity predictions based on the bulk phonon properties combined with a boundary scattering model do not capture the low thickness behavior. To match the two-dimensional lattice dynamics and molecular dynamics predictions for larger thicknesses, the three-dimensional lattice dynamics calculations require a finite specularity parameter that in some cases approaches unity. These findings point to the challenges associated with interpreting experimental thermal conductivity measurements of ultrathin silicon films, where surface roughness and a native oxide layer impact phonon transport.

Mapping phonon modes from reduced-dimensional to bulk systems

A lattice-dynamics-based algorithm for mapping the true phonon modes of a film, defined by a two-dimensional (2D) Brillouin zone, to the modes of the corresponding bulk material, defined by a three-dimensional (3D) Brillouin zone, is proposed. The mapping allows for an assessment of the common assumption that film phonon modes have the same harmonic properties as the corresponding bulk phonon modes. The algorithm is based on normal mode decomposition and is inspired by the observation that the atomic trajectories generated by the 2D eigenvectors lead to standing-wave-like behaviors in the cross-plane direction. This behavior enables an unfolding scheme that does not require an assumption of cross-plane periodicity, as used in previous approaches for mapping phonon modes. The algorithm is applied to films between two and ten unit cells thick built from Lennard-Jones (LJ) argon, whose bulk is isotropic, and graphene, whose bulk (graphite) is anisotropic. For LJ argon, the density of states deviates from that of the bulk as the film gets thinner due to phonon frequencies that shift to lower values. This shift is a result of transverse branch splitting due to the film’s anisotropy and the emergence of a quadratic acoustic branch. As such, while the mapping algorithm works well for the thicker LJ argon films, it does not perform as well for the thinner films as there is a weaker correspondence between the 2D and 3D modes. For graphene, the density of states of even the thinnest films closely matches that of graphite due to the inherent anisotropy, except for a small shift at low frequencies. As a result, the mapping algorithm works well for all thicknesses of the graphene films, indicating a strong correspondence between the 2D and 3D modes.

Non-empirical weighted Langevin mechanics for the potential escape problem: Parallel algorithm and application to the Argon clusters

Recently a non-empirical stochastic walker algorithm has been developed to search for the minimum-energy escape paths (MEP) from the minima of the potential surface (Akashi and Nagornov, 2018). This algorithm is novel in that it tracks the MEP monotonically and does not use the whole Hessian matrix but only gradient and Laplacian of the potential. In this work, we implement an parallelized version of this algorithm in a simple way. We also explore efficient ways to reduce the number of walkers required for the accurate tracking of the MEP and generate initial positions automatically. We apply the whole scheme to the Lennard-Jones argon cluster with 7–38 atoms to demonstrate the successful tracking of the reaction paths. This achievement paves the path to non-empirical simulation of rare reactions without coarse-graining or artificial potential.

The pressure in interfaces having cylindrical geometry.

While much work has been reported on the statistical mechanics and molecular simulation of interfaces of planar and spherical geometries, very little has been published on the interfaces of cylindrical geometry. The cylindrical geometry is important for the study of cylindrical micelles and particularly for nano-phases confined within cylindrical pores since the most well-defined porous materials (e.g., carbon and silicon nanotubes, SBA-15 and KIT-6 silicas) that are presently available are of this geometry. In this work, we derive the statistical mechanical equations for the pressure tensor for an interfacial region of cylindrical geometry via the virial route and for the condition of mechanical (hydrostatic) equilibrium. We also report the equation for the surface tension via the mechanical route. Monte Carlo and molecular dynamics simulation results are obtained for two example systems involving a fluid nano-phase of Lennard-Jones argon: a gas-liquid interface of cylindrical geometry and a confined nano-phase within a cylindrical carbon pore. All three diagonal elements of the pressure tensor are reported in each case, the component normal to the interface, PN = Pρρ , and the two tangential components PTϕ = Pϕϕ and PTz = Pzz , where (ρ, z, ϕ) are the usual cylindrical polar coordinates. For the cylindrical pore, the tangential pressures, Pϕϕ and Pzz , show strong compression in the adsorbed layers, as has been found in slit-shaped and spherical pores.

Quantum mechanical prediction of four-phonon scattering rates and reduced thermal conductivity of solids

Recently, first principle-based prediction of lattice thermal conductivity $\kappa$ from the perturbation theory has achieved significant success. However, it only includes three-phonon scattering due to the assumption that four-phonon and higher-order processes are generally unimportant. Also, directly evaluating the scattering rates of four-phonon and higher-order processes has been a long-standing challenge. In this work, however, we have developed a formalism to explicitly determine quantum mechanical scattering probability matrices for four-phonon scattering in the full Brillouin Zone, and by mitigating the computational challenge we have directly calculated four-phonon scattering rates. We find that four-phonon scattering rates are comparable to three-phonon scattering rates at medium and high temperatures, and they increase quadratically with temperature. As a consequence, $\kappa$ of Lennard-Jones argon is reduced by more than 60% at 80 K when four-phonon scattering is included. Also, in less anharmonic materials -- diamond, silicon, and germanium, $\kappa$ is still reduced considerably at high temperature by four-phonon scattering. Also, the thermal conductivity of optical phonons is dominated by the fourth and higher orders phonon scattering even at low temperature.

Systematic coarse-graining in nucleation theory.

In this work, we show that the standard method to obtain nucleation rate-predictions with the aid of atomistic Monte Carlo simulations leads to nucleation rate predictions that deviate 3 - 5 orders of magnitude from the recent brute-force molecular dynamics simulations [Diemand et al., J. Chem. Phys. 139, 074309 (2013)] conducted in the experimental accessible supersaturation regime for Lennard-Jones argon. We argue that this is due to the truncated state space the literature mostly relies on, where the number of atoms in a nucleus is considered the only relevant order parameter. We here formulate the nonequilibrium statistical mechanics of nucleation in an extended state space, where the internal energy and momentum of the nuclei are additionally incorporated. We show that the extended model explains the lack in agreement between the molecular dynamics simulations by Diemand et al. and the truncated state space. We demonstrate additional benefits of using the extended state space; in particular, the definition of a nucleus temperature arises very naturally and can be shown without further approximation to obey the fluctuation law of McGraw and LaViolette. In addition, we illustrate that our theory conveniently allows to extend existing theories to richer sets of order parameters.

An EQT-cDFT approach to determine thermodynamic properties of confined fluids

We present a continuum-based approach to predict the structure and thermodynamic properties of confined fluids at multiple length-scales, ranging from a few angstroms to macro-meters. The continuum approach is based on the empirical potential-based quasi-continuum theory (EQT) and classical density functional theory (cDFT). EQT is a simple and fast approach to predict inhomogeneous density and potential profiles of confined fluids. We use EQT potentials to construct a grand potential functional for cDFT. The EQT-cDFT-based grand potential can be used to predict various thermodynamic properties of confined fluids. In this work, we demonstrate the EQT-cDFT approach by simulating Lennard-Jones fluids, namely, methane and argon, confined inside slit-like channels of graphene. We show that the EQT-cDFT can accurately predict the structure and thermodynamic properties, such as density profiles, adsorption, local pressure tensor, surface tension, and solvation force, of confined fluids as compared to the molecular dynamics simulation results.

Origins of thermal conductivity changes in strained crystals

The strain-dependent phonon properties and thermal conductivities of a soft system [Lennard-Jones (LJ) argon] and a stiff system (silicon modeled using first-principles calculations) are predicted using lattice dynamics calculations and the Boltzmann transport equation. As is commonly assumed for materials under isotropic strain, the thermal conductivity of LJ argon decreases monotonically as the system moves from compression into tension. The reduction in thermal conductivity is attributed to decreases in both the phonon lifetimes and group velocities. The thermal conductivity of silicon, however, is constant in compression and only begins to decrease once the system is put in tension. The silicon lifetimes show an anomalous behavior, whereby they increase as the system moves from compression into tension, which is explained by examining the potential energy surface felt by an atom. The results emphasize the need to separately consider the harmonic and anharmonic effects of strain on material stiffness, phonon properties, and thermal conductivity.

Calculating the effect of vacancies on the thermodynamic functions and frequencies of normal vibrations of a Lennard-Jones monoatomic crystal

Four main functions of a Lennard-Jones defect argon crystal with an FCC lattice are considered in the context of the lattice gas model: Helmholtz free energy, entropy, internal energy, and heat capacity at a constant volume (allowing for normal vibrations of a solid). Properties of the defect crystal are calculated from the distribution function of the frequencies of an ideal crystal, and corrections to it that reflect local atomic vibrations around vacancies, in the context of the Lifshits linear approximation according to vacancy density. To find the free energy of a defect-free ideal crystal, frequencies of normal vibrations of the crystal are calculated with allowance for the interactions of the four nearest neighbors. The nonlocality of the chemical potentials of the atoms of a solid is discussed, and the influence of contributions from various segments of the vibrational spectrum to the values of thermodynamic functions is investigated. It is shown that ignoring the acoustic or antiphase segments of the spectrum when calculating the free energy leads to increasing of its deviation from the one calculated using the full vibrational spectrum with an increase in temperature. It is concluded that the nonequilibrium state of the defect crystal can lead to negative values of heat capacity at a constant volume.

Comparison and Evaluation of Spectral Energy Methods for Predicting Phonon Properties

Two frequency-domain methods for predicting phonon frequencies and lifetimes using the phonon spectral energy density are described. Both methods draw input from molecular dynamics simulations and lattice dynamics calculations, but differ in the form of the phonon spectral energy density. One phonon spectral energy density expression (referred to as Φ) can be formally derived from lattice dynamics theory. A similar approach in the time domain has been validated [Turney et al. Phys. Rev. B 79, 224305 (2009)]. The other phonon spectral energy density expression (referred to as Φ′) has been proposed [Thomas et al., Phys. Rev. B 81, 081411(R) (2010)] but not validated. The expressions for Φ and Φ′ are presented and then applied to predict the phonon properties and thermal conductivities of three systems: LennardJones argon, Stillinger-Weber silicon, and a carbon nanotube modeled using the Reactive Empirical Bond Order potential. Φ′ does not capture the total phonon spectral energy density predicted by Φ and therefore cannot correctly predict the phonon lifetimes or thermal conductivity. Its use in future work is discouraged and we recommend the use of Φ.

Predicting alloy vibrational mode properties using lattice dynamics calculations, molecular dynamics simulations, and the virtual crystal approximation

The virtual crystal (VC) approximation for mass disorder is evaluated by examining two model alloy systems: Lennard-Jones argon and Stillinger-Weber silicon. In both material systems, the perfect crystal is alloyed with a heavier mass species up to equal concentration. The analysis is performed using molecular dynamics simulations and lattice dynamics calculations. Mode frequencies and lifetimes are first calculated by treating the disorder explicitly and under the VC approximation, with differences found in the high-concentration alloys at high frequencies. Notably, the lifetimes of high-frequency modes are underpredicted using the VC approximation, a result we attribute to the neglect of higher-order terms in the model used to include point-defect scattering. The mode properties are then used to predict thermal conductivity under the VC approximation. For the Lennard-Jones alloys, where high-frequency modes make a significant contribution to thermal conductivity, the high-frequency lifetime underprediction leads to an underprediction of thermal conductivity compared to predictions from the Green-Kubo method, where no assumptions about the thermal transport are required. Based on observations of a minimum mode diffusivity, we propose a correction that brings the VC approximation thermal conductivities into better agreement with the Green-Kubo values. For the Stillinger-Weber alloys, where the thermal conductivity is dominated by low-frequency modes, the high-frequency lifetime underprediction does not affect the thermal conductivity prediction and reasonable agreement is found with the Green-Kubo values.

A combined quasi-continuum/Langevin equation approach to study the self-diffusion dynamics of confined fluids

In this work, we combine our earlier proposed empirical potential based quasi-continuum theory, (EQT) [A. V. Raghunathan, J. H. Park, and N. R. Aluru, J. Chem. Phys. 127, 174701 (2007)], which is a coarse-grained multiscale framework to predict the static structure of confined fluids, with a phenomenological Langevin equation to simulate the dynamics of confined fluids in thermal equilibrium. An attractive feature of this approach is that all the input parameters to the Langevin equation (mean force profile of the confined fluid and the static friction coefficient) can be determined using the outputs of the EQT and the self-diffusivity data of the corresponding bulk fluid. The potential of mean force profile, which is a direct output from EQT is used to compute the mean force profile of the confined fluid. The density profile, which is also a direct output from EQT, along with the self-diffusivity data of the bulk fluid is used to determine the static friction coefficient of the confined fluid. We use this approach to compute the mean square displacement and survival probabilities of some important fluids such as carbon-dioxide, water, and Lennard-Jones argon confined inside slit pores. The predictions from the model are compared with those obtained using molecular dynamics simulations. This approach of combining EQT with a phenomenological Langevin equation provides a mathematically simple and computationally efficient means to study the impact of structural inhomogeneity on the self-diffusion dynamics of confined fluids.

A general method for spatially coarse-graining Metropolis Monte Carlo simulations onto a lattice

A recently introduced method for coarse-graining standard continuous Metropolis Monte Carlo simulations of atomic or molecular fluids onto a rigid lattice of variable scale [X. Liu, W. D. Seider, and T. Sinno, Phys. Rev. E 86, 026708 (2012)] is further analyzed and extended. The coarse-grained Metropolis Monte Carlo technique is demonstrated to be highly consistent with the underlying full-resolution problem using a series of detailed comparisons, including vapor-liquid equilibrium phase envelopes and spatial density distributions for the Lennard-Jones argon and simple point charge water models. In addition, the principal computational bottleneck associated with computing a coarse-grained interaction function for evolving particle positions on the discretized domain is addressed by the introduction of new closure approximations. In particular, it is shown that the coarse-grained potential, which is generally a function of temperature and coarse-graining level, can be computed at multiple temperatures and scales using a single set of free energy calculations. The computational performance of the method relative to standard Monte Carlo simulation is also discussed.

Semi-infinite boundary conditions for the simulation of interfaces: The Ar/CO2(s) model revisited

We propose a method to account for the long tail corrections of dispersive forces in inhomogeneous systems. This method deals separately with the two interfaces that are usually present in a simulation setup, effectively establishing semi-infinite boundary conditions that are appropriate for the study of the interface between two infinite bulk phases. Using the wandering interface method, we calculate surface free energies of vapor-liquid, wall-liquid, and wall-vapor interfaces for a model of Lennard-Jones argon adsorbed on solid carbon dioxide. The results are employed as input to Young's equation, and the wetting temperature located. This estimate is compared with predictions from the method of effective interface potentials and good agreement is found. Our results show that truncating Ar-Ar interactions at two and a half molecular diameters results in a dramatic decrease of the wetting temperature of about 40%.

Anderson Localization of Phonons in Random Multilayer Thin Films

ABSTRACT1020 Joules of energy are generated by the United States each year; 60% of this energy is lost to waste heat [1]. Thermoelectric based energy scavenging has tremendous potential for the recovery of significant quantities of this waste heat. However, utilization of thermoelectric devices is limited due to relatively low energy conversion efficiency and the utilization of relatively scarce materials. This work focuses on generating sustainable and efficient thermoelectric materials through modifications to the lattice vibrations of materials with excellent thermoelectric electronic properties (Seebeck coefficients larger than 500 μV/K). In particular, Anderson localization of phonons in random multilayer thin films has been explored as a means for reducing lattice thermal conductivity to values approaching that of aerogels (∼10 mW/m-K). Silicon has been a sample of choice due to its high crust abundance and Seebeck coefficient. Reverse non-equilibrium molecular dynamics simulations have been utilized to determine the thermal conductivity of structures of interest. Simulations with pure Lennard-Jones argon solids have been performed to establish a methodology and to characterize the effect of different kinds of disorder prior to the examination of silicon. The simulation results indicate that mass disorder confined to randomly selected planes to be an effective way in which to reduce lattice thermal conductivity with the lattice thermal conductivity decreasing by a factor of thirty (to 4 mW/m-K) in the argon case and a factor of over ten thousand (to 15 mW/m-K) for silicon. Based on models in which the charge carrier mean free path is limited by scattering from the planes with mass disorder, the mobility of silicon is expected to reach values of 10 cm2/V-s. At this mobility the thermoelectric figure of merit, ZT, (utilizing the Wiedeman-Franz law to calculate the electronic thermal conductivity) varies between 4.5 and 11 as the mass ratio of the disordered planes is varied from 4 to 10 in 20% of the lattice planes. These results indicate that the pursuit of nanostructured thermoelectric materials in the form of random multilayers may provide a path to efficient and sustainable thermoelectric materials.

Size effects in molecular dynamics thermal conductivity predictions

We predict the bulk thermal conductivity of Lennard-Jones argon and Stillinger-Weber silicon using the Green-Kubo (GK) and direct methods in classical molecular dynamics simulations. While system-size-independent thermal conductivities can be obtained with less than 1000 atoms for both materials using the GK method, the linear extrapolation procedure [Schelling et al., Phys. Rev. B 65, 144306 (2002)] must be applied to direct method results for multiple system sizes. We find that applying the linear extrapolation procedure in a manner consistent with previous researchers can lead to an underprediction of the GK thermal conductivity (e.g., by a factor of 2.5 for Stillinger-Weber silicon at a temperature of 500 K). To understand this discrepancy, we perform lattice dynamics calculations to predict phonon properties and from these, length-dependent thermal conductivities. From these results, we find that the linear extrapolation procedure is only accurate when the minimum system size used in the direct method simulations is comparable to the largest mean-free paths of the phonons that dominate the thermal transport. This condition has not typically been satisfied in previous works. To aid in future studies, we present a simple metric for determining if the system sizes used in direct method simulations are sufficiently large so that the linear extrapolation procedure can accurately predict the bulk thermal conductivity.