Nonequilibrium Molecular Dynamics Algorithm Research Articles

A Reverse Nonequilibrium Molecular Dynamics Algorithm for Coupled Mass and Heat Transport in Mixtures.

We present a new method for introducing stable nonequilibrium concentration gradients in molecular dynamics simulations of mixtures. This method extends earlier reverse nonequilibrium molecular dynamics (RNEMD) methods, which use kinetic energy scaling moves to create temperature or velocity gradients. In the new scaled particle flux (SPF-RNEMD) algorithm, energies and forces are computed simultaneously for a molecule existing in two nonadjacent regions of a simulation box, and the system evolves under a linear combination of these interactions. A continuously increasing particle scaling variable is responsible for the migration of the molecule between the regions as the simulation progresses, allowing for simulations under an applied particle flux. To test the method, we investigate diffusivity in mixtures of identical but distinguishable particles and in a simple mixture of multiple Lennard-Jones particles. The resulting concentration gradients provide Fick diffusion constants for mixtures. We also discuss using the new method to obtain coupled transport properties using simultaneous particle and thermal fluxes to compute the temperature dependence of the diffusion coefficient and activation energies for diffusion from a single simulation. Lastly, we demonstrate the use of this new method in interfacial systems by computing the diffusive permeability of a molecular fluid moving through a nanoporous graphene membrane.

A low‐disturbance nonequilibrium molecular dynamics algorithm applied to the determination of thermal conductivities

A new nonequilibrium molecular dynamics algorithm is proposed for the determination of thermal conductivity and other transport properties. The proposed algorithm aims at diminishing the energy drift problem observed in this type of method while conserving linear momentum and being compatible with constrained molecules. The features of the proposed algorithm are evaluated by determining thermal conductivities of water at 323 K, n‐octane at 300 K, and argon close to its triple point, and by comparing these results with those obtained using established methods. The analysis of systems presenting diverse molecular characteristics allowed us to assess the usefulness of the proposed algorithm. The energy drift and temperature variation were reduced in the range of 10–80%, depending on the parameters of the proposed algorithm and the characteristics of the system. The determined thermal conductivities showed good agreement when compared to experimental and simulation data. © 2015 American Institute of Chemical Engineers AIChE J, 61: 2881–2890, 2015

A reverse nonequilibrium molecular dynamics method for calculating the mutual diffusion coefficient for binary fluids

This work introduces a nonequilibrium molecular dynamics method for calculating the mutual diffusion coefficient for mixtures. The method is based on the idea of the reverse nonequilibrium molecular dynamics algorithms, and artificially generates a mass flux through the Lennard–Jones mixtures by suitably exchanging particle positions and velocities in different regions. The analysis of the resulting steady-state concentration profiles allows the calculation of the mutual diffusion coefficient. As a test, this method is applied to the calculation of the mutual diffusion coefficient of Ar/Kr systems, and reasonable results are obtained. At the same time, it is observed that the mutual diffusion coefficient strongly depends on the temperature and composition of the mixture. The method can be easily extended to other fluid mixtures and be adopted when studying the nonequilibrium fluid mixtures.

Thermal conductivity calculation of complex (dusty) plasmas

The thermal conductivity of three-dimensional (3D) strongly coupled complex (dusty) plasmas has been calculated through the improved Evan-Gillan nonequilibrium molecular dynamics (NEMD) algorithm. The extensive NEMD simulations are performed to study the performance of the algorithm and compared the results determined for perturbed heat energy current to the results obtained by equilibrium molecular dynamics (EMD) simulations. The calculations show that the present algorithm gives accurate results with fast convergence and small size effects over a wide range of plasma coupling and screening parameters. The present simulation results are in agreement with part of others NEMD and EMD data in the literature with simulation values generally overpredicting the thermal conductivity by 3%–20%, depending on plasma parameters. It is shown that the homogenous perturbed method can be employed to estimate the thermal conductivity and to understand the fundamental behaviors in 3D complex Yukawa liquids.

Nanoscale simple-fluid behavior under steady shear

In this study, we use two nonequilibrium molecular dynamics algorithms, boundary-driven shear and homogeneous shear, to explore the rheology and flow properties of a simple fluid undergoing steady simple shear. The two distinct algorithms are designed to elucidate the influences of nanoscale confinement. The results of rheological material functions, i.e., viscosity and normal pressure differences, show consistent Newtonian behaviors at low shear rates from both systems. The comparison validates that confinements of the order of 10 nm are not strong enough to deviate the simple fluid behaviors from the continuum hydrodynamics. The non-Newtonian phenomena of the simple fluid are further investigated by the homogeneous shear simulations with much higher shear rates. We observe the "string phase" at high shear rates by applying both profile-biased and profile-unbiased thermostats. Contrary to other findings where the string phase is found to be an artifact of the thermostats, we perform a thorough analysis of the fluid microstructures formed due to shear, which shows that it is possible to have a string phase and second shear thinning for dense simple fluids.

Twist viscosities and flow alignment of biaxial nematic liquid crystal phases of a soft ellipsoid-string fluid studied by molecular dynamics simulation

We have calculated the twist viscosity and the alignment angle between the director and the stream lines in shear flow of a liquid crystal model system, which forms biaxial nematic liquid crystals, as functions of the density, from the Green-Kubo relations by equilibrium molecular dynamics simulation and by a nonequilibrium molecular dynamics algorithm, where a torque conjugate to the director angular velocity is applied to rotate the director. The model system consists of a soft ellipsoid-string fluid where the ellipsoids interact according a repulsive version of the Gay-Berne potential. Four different length-to-width-to-breadth ratios have been studied. On compression, this system forms discotic or calamitic uniaxial nematic phases depending on the dimensions of the molecules, and on further compression a biaxial nematic phase is formed. In the uniaxial nematic phase there is one twist viscosity and one alignment angle. In the biaxial nematic phase there are three twist viscosities and three alignment angles corresponding to the rotation around the various directors and the different alignments of the directors relative to the stream lines, respectively. It is found that the smallest twist viscosity arises by rotation around the director formed by the long axes, the second smallest one arises by rotation around the director formed by the normals of the broadsides, and the largest one by rotation around the remaining director. The first twist viscosity is rather independent of the density whereas the last two ones increase strongly with density. One finds that there is one stable director alignment relative to the streamlines, namely where the director formed by the long axes is almost parallel to the stream lines and where the director formed by the normals of the broadsides is almost parallel to the shear plane. The relative magnitudes of the components of the twist viscosities span a fairly wide interval so this model should be useful for parameterisation experimental data.

Planar mixed flow and chaos: Lyapunov exponents and the conjugate-pairing rule

In this work we characterize the chaotic properties of atomic fluids subjected to planar mixed flow, which is a linear combination of planar shear and elongational flows, in a constant temperature thermodynamic ensemble. With the use of a recently developed nonequilibrium molecular dynamics algorithm, compatible and reproducible periodic boundary conditions are realized so that Lyapunov spectra analysis can be carried out for the first time. Previous studies on planar shear and elongational flows have shown that Lyapunov spectra organize in different ways, depending on the character of the defining equations of the system. Interestingly, planar mixed flow gives rise to chaotic spectra that, on one hand, contain elements common to those of shear and elongational flows but also show peculiar, unique traits. In particular, the influence of the constituent flows in regards to the conjugate-pairing rule (CPR) is analyzed. CPR is observed in homogeneously thermostated systems whose adiabatic (or unthermostated) equations of motion are symplectic. We show that the component associated with the shear tends to selectively excite some of those degrees, and is responsible for violations in the rule.

Examining different NEMD methods in simulating nanoscale fluid at high shear rates

In this paper, the non-Newtonian behaviour of a simple fluid is simulated via three different non-equilibrium molecular dynamics (NEMD) algorithms: boundary-driven shear; homogeneous shear; and reverse-NEMD (RNEMD). The pros and cons of each approach are discussed and results are compared. It is found that the accessible shear rate is very limited when using boundary-driven shear and RNEMD, whereas homogeneous shear is capable of imposing a substantially higher shear rate. However, the boundary-driven shear reproduces rheological experiment setups and the RNEMD yields much better statistically calculated properties. The material functions, velocity profiles, and microstructures of the fluid are examined and the formations of fluid structures are elucidated. The homogeneous shear is shown to be the most appropriate NEMD method to simulate fluids under strong shear, with the smallest number of limitations.

Nonequilibrium molecular dynamics calculation of the thermal conductivity based on an improved relaxation scheme

A nonequilibrium molecular dynamics (NEMD) method using stochastic energy injection and removal as uniform heat sources and sinks is developed to calculate the thermal conductivity. The stochastic energy is generated by a Maxwell function generator and is imposed on only a few individual molecules each time step. The relaxation of the thermal perturbation is improved compared to other NEMD algorithms because there are no localized heat source and sink slab regions in the system. The heat sources are uniformly distributed in the right half of the system while the sinks are in the left half, which leads to a periodically quadratic temperature distribution that is almost sinusoidal. The thermal conductivity is then easily calculated from the mean temperatures of the right and left half systems rather than by fitting the temperature profiles. This improved relaxation NEMD scheme is used to calculate the thermal conductivities of liquid and solid argons. It shows that the present algorithm gives accurate results with fast convergence and small size effects. Other stochastic energy perturbation, e.g., thermal noise, can be used to replace the Maxwell-type perturbation used in this paper to make the improved relaxation scheme more effective.

Equilibrium and nonequilibrium molecular dynamics simulations of the thermal conductivity of molten alkali halides

The thermal conductivity of molten NaCl and KCl was calculated through the Evans-Gillan nonequilibrium molecular dynamics (NEMD) algorithm and Green-Kubo equilibrium molecular dynamics (EMD) simulations. The EMD simulations were performed for a "binary" ionic mixture and the NEMD simulations assumed a pure system for reasons discussed in this work. The cross thermoelectric coefficient obtained from Green-Kubo EMD simulations is discussed in terms of the homogeneous thermoelectric power or Seebeck coefficient of these materials. The thermal conductivity obtained from NEMD simulations is found to be in very good agreement with that obtained through Green-Kubo EMD simulations for a binary ionic mixture. This result points to a possible cancellation between the neglected "partial enthalpy" contribution to the heat flux associated with the interdiffusion of one species through the other and that part of the thermal conductivity related to the coupled fluxes of charge and heat in "binary" ionic mixtures.

Homogeneous non-equilibrium molecular dynamics simulations of viscous flow: techniques and applications

We provide a review of the literature for non-equilibrium molecular dynamics (NEMD) simulations of homogeneous fluids. Our review focuses on techniques for simulations of shear and elongational flows in viscous fluids and covers the formulation and application of NEMD algorithms for atomic and molecular fluids. We provide a set of expositions that can be effectively used as guidelines to formulate the relevant equations of motion, periodic boundary conditions and thermostats. We also provide a survey of applications in a convenient tabular form as an aid to researchers who wish to use NEMD to study transport phenomena.

Rheological and structural studies of liquid decane, hexadecane, and tetracosane under planar elongational flow using nonequilibrium molecular-dynamics simulations

We report for the first time rheological and structural properties of liquid decane, hexadecane, and tetracosane using nonequilibrium molecular-dynamics (NEMD) simulations under planar elongational flow (PEF). The underlying NEMD algorithm employed is the so-called p-SLLOD algorithm [C. Baig, B. J. Edwards, D. J. Keffer, and H. D. Cochran, J. Chem. Phys. 122, 114103 (2005)]. Two elongational viscosities are measured, and they are shown not to be equal to each other, indicating two independent viscometric functions in PEF. With an appropriate definition, it is observed that the two elongational viscosities converge to each other at very low elongation rates, i.e., in the Newtonian regime. For all three alkanes, tension-thinning behavior is observed. At high elongation rates, chains appear to be fully stretched. This is supported by the result of the mean-square end-to-end distance of chains (R(ete2)) and the mean-square radius of gyration of chains (R(g2)), and further supported by the result of the intramolecular Lennard-Jones (LJ) potential energy. It is also observed that (R(ete2)) and (R(g2)) show a different trend as a function of strain rate from those in shear flow: after reaching a plateau value, (R(ete2)) and (R(g2)) are found to increase further as elongation rate increases. A minimum in the hydrostatic pressure is observed for hexadecane and tetracosane at about epsilon(msigma2/epsilon)1/2=0.02. This phenomenon is shown to be associated with the intermolecular LJ potential energy. The bond-bending and torsional energies display similar trends, but a different behavior is observed for the bond-stretching energy. An important observation common in these three bonded-intramolecular interactions is that all three modes are suppressed to a small value at high elongation rates. We conjecture that a liquid-crystal-like, nematic structure is present in these systems at high elongation rates, which is characterized by a strong chain alignment with a fully stretched conformation.

Computing the Soret coefficient in aqueous mixtures using boundary driven nonequilibrium molecular dynamics

We have computed the Soret coefficient in aqueous mixtures using a boundary driven nonequilibrium molecular dynamics algorithm and standard molecular force fields. The choice of this specific approach is justified by the nature of the mixtures studied here. Four aqueous solutions, including methanol, ethanol, acetone, and dimethyl-sulfoxide (DMSO) have been studied at ambient conditions for different compositions. The experimental behavior of water-alcohol mixtures was reproduced, including the change of sign of the Soret coefficient with composition, in excellent agreement with existing experimental data. The methodology has been applied to obtain pure predictions for water-acetone and water-DMSO where no experimental data are accessible. A change of sign is also observed in the same range of composition as in water-alcohol mixtures. It is suggested that the nature and strength of the molecular interactions, rather than the mass or shape ratio of the components, dominates the behavior of the Soret coefficient versus composition for the aqueous associating mixtures studied here.

A nonequilibrium molecular dynamics method for thermal conductivities based on thermal noise

We developed a nonequilibrium molecular dynamics (NEMD) method for calculating thermal conductivities. In contrast to other NEMD algorithms, here only the heat sink is localized, whereas the heat source can be uniformly distributed throughout the system. The noise due to cutting off the pair forces or to integration errors is such a uniform heat source. In traditional NEMD methods it is normally considered a nuisance factor. The new algorithm accounts for it and uses it. The algorithm is easy to derive, analyse and implement. Moreover, it circumvents the need to calculate energy fluxes. It is tested on the enhanced simple-point charge model for liquid water and reproduces the known thermal conductivity of this model liquid of 0.81 W m(-1) K(-1). It can be generalized to situations, where the thermal noise is replaced by another uniform heat source, or to the inverse situation, where the heat source is localized but the heat sink extends over the entire system.

Molecular Dynamics Simulations of Multicomponent Diffusion. 2. Nonequilibrium Method

We describe a new nonequilibrium molecular dynamics algorithm for simulating multicomponent diffusion within the Stefan−Maxwell framework. The approach developed here allows for the efficient computation of a full diffusion-coefficient matrix for solutions containing any number of species. Due to our interest in electrolytes, the method is tested by simulating aqueous KCl and NaCl salt solutions in the concentration range 1−4 m. The results are in semiquantitative agreement with experiment and in quantitative agreement with previous Green−Kubo calculations for our potential models. This provides strong evidence that the equilibrium and nonequilibrium methods are correct.

Thermal diffusion in Lennard–Jones fluids in the frame of the law of the corresponding states

This work is related to the definition of a reduced thermal diffusion coefficient thanks to numerical microscale molecular dynamics simulations. This cross transport process, also called Soret effect, couples mass flux and thermal gradient and is still largely misunderstood. For this study, we have applied a boundary driven non-equilibrium molecular dynamics algorithm on Lennard–Jones spheres mixtures. Simulations have been performed at a constant reduced supercritical state, using a van der Waals’ one fluid approximation in order to fulfil the law of the corresponding states. In binary mixtures, we have studied the molecular parameters and the molar fraction influences on thermal diffusion separately and then combined. It is shown that, on pressure and on thermal conductivity, the corresponding states law is fulfilled for a wide range of molecular parameters ratios. In this frame, we have then constructed simple correlations which relate thermal diffusion factor to the mixture parameters. Combining the relations obtained, a reduced thermal diffusion factor taking into account all the various contributions has been defined. Finally, it is shown that this relation enables us to estimate thermal diffusion in various binary and ternary mixtures of Lennard–Jones spheres representing alkanes with a maximum deviation of 15%.

Shear Viscosity of Benzene, Toluene, and p-Xylene by Non-equilibrium Molecular Dynamics Simulations

Green and Kubo 1 showed that the phenomenological coefficients describing many transport processes and timedependent phenomena in general could be written as integrals over a certain type of function called a timecorrelation function. The Green-Kubo formulas are the formal expressions for hydrodynamic field variables and some of the thermodynamic properties in terms of the microscopic variables of an N-particle system. The identification of microscopic expressions for macroscopic variables is made by a process of comparison of the conservation equations of hydrodynamics with the microscopic equations of change for conserved densities. The importance of these formulas is three-fold: they provide an obvious method for calculating transport coefficients using computer simulation, a convenient starting point for constructing analytic theories for non-equilibrium processes, and an essential information for designing non-equilibrium molecular dynamics (NEMD) algorithm. The Green-Kubo formula for the shear viscosity is given by

Thermal diffusion sensitivity to the molecular parameters of a binary equimolar mixture, a non-equilibrium molecular dynamics approach

The goal of this study is to analyse how the thermal diffusion process is dependent on molecular parameters when describing a fluid mixture. To estimate the associated transport coefficient, which is the thermal diffusion factor α T, a non-equilibrium molecular dynamics algorithm has been applied on equimolar binary mixtures of Lennard–Jones (LJ) particles in supercritical conditions. Firstly, it is shown that this model is able to correctly estimate α T for simple alkane mixtures, provided there are a sufficient number of particles and long enough simulations. Then, using various mixing rules, the separate influences of the mass, the moment of inertia, the atomic diameter and the interaction strength have been studied. Results indicate that the molar fraction of the component, having the smallest mass and moment of inertia as well as the biggest radius and the strongest potential, tends to increase in the hot area. Elsewhere, simulations for various cross-interaction parameters show that α T is extremely sensitive to the intermolecular pair potential between unlike particles. Finally, results on methane/normal alkane mixtures indicate that a simple sum between the separate contributions provides a reliable evaluation of α T only when the molecular parameter ratios between the two components are close to 1.

On thermal diffusion in binary and ternary Lennard-Jones mixtures by non-equilibrium molecular dynamics

Molecular simulation appears to be an alternative to experiment for the estimation of transport and thermodynamics properties of fluid mixtures, which is of primary importance in the evaluation of the initial state of a petroleum reservoir. In this study, a non-equilibrium molecular dynamics algorithm has been applied to mixtures of Lennard-Jones spheres in order to compute the thermal diffusion process. The pertinence of such an approach to simple alkane mixtures is shown. The separate influences on the thermal diffusion of the molecular features in binary equimolar mixtures are then summarized. Simulations on binary non-equimolar mixtures have been performed as well. The results indicate an increase in the thermal diffusion process with increasing molar fraction of the lightest component. Moreover, this increase is enhanced with increasing difference in the number of carbons between the two alkanes. Then, a simple method, which yields results consistent with simulations, is proposed to predict thermal diffusion for the whole range of molar fractions starting only from the equimolar value. Finally, for ternary mixtures, the law of the corresponding states is shown to be valid when the appropriate mixing rules are applied, which allows the estimation of thermal diffusion in such mixtures from equivalent binary mixtures.

Correspondence between configurational temperature and molecular kinetic temperature thermostats

Molecular fluids undergoing shear flow are often modeled using a homogeneous nonequilibrium molecular dynamics algorithm. To reach a steady state, this method must be used in conjunction with a thermostating mechanism which duplicates the heat dissipation in the experimental setup (e.g., by conduction to the shearing boundaries). The most commonly used type of thermostat involves fixing the center of mass kinetic (c.m.) temperature. Though perfectly valid, this approach does not seem to be the most realistic for a molecular fluid since heat is removed only through the 3 degrees of freedom of the center of mass for each molecule. The second type of thermostat involves fixing the “atomic” kinetic temperature and therefore takes into account all degrees of freedom. However, since the streaming velocity of atoms within their constituent molecules is unknown, the implementation of such a thermostat is problematic and relies on incorrect assumptions on the streaming velocity of atoms. The recently developed configurational temperature thermostat requires no assumption on the streaming velocity of atoms and takes into account all degrees of freedom. Using a configurational temperature thermostat to thermostat homogeneous shear flow thus seems to be a more realistic approach than the c.m. kinetic thermostat. In this work, we apply this configurational temperature thermostat to the study of linear alkanes (C10 and C20) undergoing shear flow. The results so obtained are compared with those obtained using a c.m. kinetic thermostat. Our aims are (1) to test the influence of the total number of degrees of freedom of the system, (2) to make a connection between the results obtained with the two types of thermostats. By carefully examining the energies of the internal modes, we have been able to characterize the loss of accuracy of a c.m. kinetic thermostat at high shear rates and for high molecular weight compounds. Finally, we establish a correspondence between the two types of thermostats by showing that, for the internal modes, a simulation at a fixed c.m. kinetic temperature is equivalent to a simulation at a fixed but higher configurational temperature.