A gentler approach to RNEMD: nonisotropic velocity scaling for computing thermal conductivity and shear viscosity

We present a new method for introducing stable nonequilibrium velocity and temperature gradients in molecular dynamics simulations of heterogeneous systems. This method extends earlier reverse nonequilibrium molecular dynamics (RNEMD) methods which use momentum exchange swapping moves. The standard swapping moves can create nonthermal velocity distributions and are difficult to use for interfacial calculations. By using nonisotropic velocity scaling (NIVS) on the molecules in specific regions of a system, it is possible to impose momentum or thermal flux between regions of a simulation while conserving the linear momentum and total energy of the system. To test the method, we have computed the thermal conductivity of model liquid and solid systems as well as the interfacial thermal conductivity of a metal-water interface. We find that the NIVS-RNEMD improves the problematic velocity distributions that develop in other RNEMD methods.     

Slip flow in graphene nanochannels

We investigate the hydrodynamic boundary condition for simple nanofluidic systems such as argon and methane flowing in graphene nanochannels using equilibrium molecular dynamics simulations (EMD) in conjunction with our recently proposed method [J. S. Hansen, B. D. Todd, and P. J. Daivis, Phys. Rev. E 84, 016313 (2011)]. We first calculate the fluid-graphene interfacial friction coefficient, from which we can predict the slip length and the average velocity of the first fluid layer close to the wall (referred to as the slip velocity). Using direct nonequilibrium molecular dynamics simulations (NEMD) we then calculate the slip length and slip velocity from the streaming velocity profiles in Poiseuille and Couette flows. The slip lengths and slip velocities from the NEMD simulations are found to be in excellent agreement with our EMD predictions. Our EMD method therefore enables one to directly calculate this intrinsic friction coefficient between fluid and solid and the slip length for a given fluid and solid, which is otherwise tedious to calculate using direct NEMD simulations at low pressure gradients or shear rates. The advantages of the EMD method over the NEMD method to calculate the slip lengths/flow rates for nanofluidic systems are discussed, and we finally examine the dynamic behaviour of slip due to an externally applied field and shear rate.     

Molecular diffusion and slip boundary conditions at smooth surfaces with periodic and random nanoscale textures

The influence of periodic and random surface textures on the flow structure and effective slip length in Newtonian fluids is investigated by molecular dynamics (MD) simulations. We consider a situation where the typical pattern size is smaller than the channel height and the local boundary conditions at wetting and nonwetting regions are characterized by finite slip lengths. In the case of anisotropic patterns, transverse flow profiles are reported for flows over alternating stripes of different wettability when the shear flow direction is misaligned with respect to the stripe orientation. The angular dependence of the effective slip length obtained from MD simulations is in good agreement with hydrodynamic predictions provided that the stripe width is larger than several molecular diameters. We found that the longitudinal component of the slip velocity along the shear flow direction is proportional to the interfacial diffusion coefficient of fluid monomers in that direction at equilibrium. In case of random textures, the effective slip length and the diffusion coefficient of fluid monomers in the first layer near the heterogeneous surface depend sensitively on the total area of wetting regions.     

Slip length of water on graphene: limitations of non-equilibrium molecular dynamics simulations

Data for the flow rate of water in carbon nanopores is widely scattered, both in experiments and simulations. In this work, we aim at precisely quantifying the characteristic large slip length and flow rate of water flowing in a planar graphene nanochannel. First, we quantify the slip length using the intrinsic interfacial friction coefficient between water and graphene, which is found from equilibrium molecular dynamics (EMD) simulations. We then calculate the flow rate and the slip length from the streaming velocity profiles obtained using non-equilibrium molecular dynamics (NEMD) simulations and compare with the predictions from the EMD simulations. The slip length calculated from NEMD simulations is found to be extremely sensitive to the curvature of the velocity profile and it possesses large statistical errors. We therefore pose the question: Can a micrometer range slip length be reliably determined using velocity profiles obtained from NEMD simulations? Our answer is "not practical, if not impossible" based on the analysis given as the results. In the case of high slip systems such as water in carbon nanochannels, the EMD method results are more reliable, accurate, and computationally more efficient compared to the direct NEMD method for predicting the nanofluidic flow rate and hydrodynamic boundary condition.     

Thermal relaxation mechanism and role of chemical functionalization in fullerene solutions

Using molecular-dynamics simulations we investigate thermal relaxation of C60 and C84 molecules suspended in octane liquid. Pristine fullerenes exhibit relatively slow relaxation due to weak thermal coupling with the liquid. A comparison of the interfacial transport characteristics obtained from relaxation simulations with those obtained from equilibrium simulations and fluctuation-dissipation theorem analysis demonstrates that the relaxation process involves two main steps: (i) energy flow from high- to low-frequency modes within the fullerene, and (ii) energy flow from low-frequency fullerene modes to the liquid. Functionalization of fullerenes with alkene chains leads to significant reduction of the thermal relaxation time. The relaxation time of functionalized fullerenes becomes independent from the functionalizing chain length beyond approximately 10 carbon segments; this can be understood in terms of thermal conductivity along the chain and heat transfer between the chain and the solvent.     

Boundary slip and wetting properties of interfaces: correlation of the contact angle with the slip length

Correlations between contact angle, a measure of the wetting of surfaces, and slip length are developed using nonequilibrium molecular dynamics for a Lennard-Jones fluid in Couette flow between graphitelike hexagonal-lattice walls. The fluid-wall interaction is varied by modulating the interfacial energy parameter epsilonr=epsilonsfepsilonff and the size parameter sigmar=sigmasfsigmaff, (s=solid, f=fluid) to achieve hydrophobicity (solvophobicity) or hydrophilicity (solvophilicity). The effects of surface chemistry, as well as the effects of temperature and shear rate on the slip length are determined. The contact angle increases from 25 degrees to 147 degrees on highly hydrophobic surfaces (as epsilonr decreases from 0.5 to 0.1), as expected. The slip length is functionally dependent on the affinity strength parameters epsilonr and sigmar: increasing logarithmically with decreasing surface energy epsilonr (i.e., more hydrophobic), while decreasing with power law with decreasing size sigmar. The mechanism for the latter is different from the energetic case. While weak wall forces (small epsilonr) produce hydrophobicity, larger sigmar smoothes out the surface roughness. Both tend to increase the slip. The slip length grows rapidly with a high shear rate, as wall velocity increases three decades from 100 to 10(5) ms. We demonstrate that fluid-solid interfaces with low epsilonr and high sigmar should be chosen to increase slip and are prime candidates for drag reduction.     

Phase field theory of interfaces and crystal nucleation in a eutectic system of fcc structure: I. Transitions in the one-phase liquid region

The phase field theory (PFT) has been applied to predict equilibrium interfacial properties and nucleation barrier in the binary eutectic system Ag-Cu using double well and interpolation functions deduced from a Ginzburg-Landau expansion that considers fcc (face centered cubic) crystal symmetries. The temperature and composition dependent free energies of the liquid and solid phases are taken from CALculation of PHAse Diagrams-type calculations. The model parameters of PFT are fixed so as to recover an interface thickness of approximately 1 nm from molecular dynamics simulations and the interfacial free energies from the experimental dihedral angles available for the pure components. A nontrivial temperature and composition dependence for the equilibrium interfacial free energy is observed. Mapping the possible nucleation pathways, we find that the Ag and Cu rich critical fluctuations compete against each other in the neighborhood of the eutectic composition. The Tolman length is positive and shows a maximum as a function of undercooling. The PFT predictions for the critical undercooling are found to be consistent with experimental results. These results support the view that heterogeneous nucleation took place in the undercooling experiments available at present. We also present calculations using the classical droplet model [classical nucleation theory (CNT)] and a phenomenological diffuse interface theory (DIT). While the predictions of the CNT with a purely entropic interfacial free energy underestimate the critical undercooling, the DIT results appear to be in a reasonable agreement with the PFT predictions.     

Calculating fifth-order Raman signals for various molecular liquids by equilibrium and nonequilibrium hybrid molecular dynamics simulation algorithms

The fifth-order two-dimensional (2D) Raman signals have been calculated from the equilibrium and nonequilibrium (finite field) molecular dynamics simulations. The equilibrium method evaluates response functions with equilibrium trajectories, while the nonequilibrium method calculates a molecular polarizability from nonequilibrium trajectories for different pulse configurations and sequences. In this paper, we introduce an efficient algorithm which hybridizes the existing two methods to avoid the time-consuming calculations of the stability matrices which are inherent in the equilibrium method. Using nonequilibrium trajectories for a single laser excitation, we are able to dramatically simplify the sampling process. With this approach, the 2D Raman signals for liquid xenon, carbon disulfide, water, acetonitrile, and formamide are calculated and discussed. Intensities of 2D Raman signals are also estimated and the peak strength of formamide is found to be only five times smaller than that of carbon disulfide.     

Interfacial properties of water + alcohol mixtures

Mixtures of water with alcohol are important in numerous engineering applications. Caused by the polarity of water and alcohol self-association of water and alcohol as well cross-association between water and alcohol appear in such complex mixtures. These features show significant impact on physical and chemical properties, especially phase equilibrium behaviour and hence interfacial properties. The Cahn–Hilliard theory was combined with original statistical associated fluid theory equation of states (SAFT EOS) in order to describe both the phase behaviour and interfacial properties with respect of association. The paper focuses on theoretical investigations of surface tension, density profiles, surface thickness in vapour–liquid or vapour–liquid–liquid equilibrium of mixtures of water with ethanol or 1-butanol. Results of vapour–liquid equilibrium surface tension calculations were compared with experimental data taken from the literature.     

Interfacial slip at the thermotropic liquid‐crystalline polymer/poly (ethylene naphthalate) interface

The unique rheological properties of a thermotropic liquid‐crystalline polymer (TLCP) were first studied. The thermal and shear history of the TLCP was found to play a critical role in its rheological properties. Crystallites were observed in the TLCP melt even above the melting temperature detected by differential scanning calorimetry. Because interfacial slip had long been suggested as an important reason for viscosity reduction in TLCP/thermoplastic blends, for the first time, interfacial slip at the TLCP/poly(ethylene naphthalate) (PEN) interface was investigated with an energy model. The model quantified the degree of interfacial slip at the TLCP/PEN interface by an energy factor. The calculated energy factors revealed a high degree of interfacial slip at the TLCP/PEN interface. It was proposed that the high rigidity of rodlike TLCP chains and their alignment parallel to the interface prevented mutual entanglements at the TLCP/PEN interface. The lack of mutual entanglements promoted the interfacial slip. © 2003 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 42: 302–315, 2004     

Static and dynamic properties of the interface between a polymer brush and a melt of identical chains

Molecular-dynamics simulations of a short-chain polymer melt between two brush-covered surfaces under shear have been performed. The end-grafted polymers which constitute the brush have the same chemical properties as the free chains in the melt and provide a soft deformable substrate. Polymer chains are described by a coarse-grained bead-spring model, which includes excluded volume and backbone connectivity of the chains. The grafting density of the brush layer offers a way of controlling the behavior of the surface without altering the molecular interactions. We perform equilibrium and nonequilibrium molecular-dynamics simulations at constant temperature and volume using the dissipative particle dynamics thermostat. The equilibrium density profiles and the behavior under shear are studied as well as the interdigitation of the melt into the brush, the orientation on different length scales (bond vectors, radius of gyration, and end-to-end vector) of free and grafted chains, and velocity profiles. The obtained boundary conditions and slip length show a rich behavior as a function of grafting density and shear velocity.     

Electrode polarization effects on interfacial kinetics of ionic liquid at graphite surface: An extended lagrangian-based constant potential molecular dynamics simulation study

Computational models including electrode polarization can be essential to study electrode/electrolyte interfacial phenomena more realistically. We present here a constant-potential classical molecular dynamics simulation method based on the extended Lagrangian formulation where the fluctuating electrode atomic charges are treated as independent dynamical variables. The method is applied to a graphite/ionic liquid system for the validation and the interfacial kinetics study. While the correct adiabatic dynamics is achieved with a sufficiently small fictitious mass of charge, static properties have been shown to be almost insensitive to the fictitious mass. As for the kinetics study, electrical double layer (EDL) relaxation and ion desorption from the electrode surface are considered. We found that the polarization slows EDL relaxation greatly whereas it has little impact on the ion desorption kinetics. The findings suggest that the polarization is essential to estimate the kinetics in nonequilibrium processes, not in equilibrium. © 2019 Wiley Periodicals, Inc.     

Single-molecule dynamics of mechanical coiled-coil unzipping

We use atomic force microscopy (AFM) to mechanically unzip and rezip a double-stranded coiled-coil structure at varying pulling velocities. We find that force-extension traces exhibit hysteresis that grows with increasing pulling velocity. This shows that coiled-coil unzipping and rezipping do not occur in thermal equilibrium on our experimental time scale. We present a nonequilibrium simulation that fully reproduces the hysteresis effects, giving detailed insight into dynamics of coiled-coil folding. Using this model, we find that seed formation is responsible for the hysteresis. The seed consists of four alpha-helical turns on both strands of the coiled coil. To obtain equilibrium information from our nonequilibrium experiments, we used the Crooks fluctuation theorem (CFT) to calculate the equilibrium free energy of folding for all of the different pulling velocities. The paper presented here lays the groundwork for the study of self-assembly properties of many physiologically relevant coiled-coil structures at the single-molecule level.     

Modeling processes of thermal decomposition of hydrocarbon fuels in heated channels

We introduce a mathematical model of the nonequilibrium process of thermal decomposition of hydrocarbon fuel in heated channels of a ramjet combustion chamber cooling system. This mathematical model is based on describing the process using intermediate asymptotics formed when taking into account the equilibrium gas composition, which is determined using open source software for calculating the equilibrium state of the chemical reaction products. A procedure was introduced allowing at different stages of the process of thermal decomposition of fuel to separate kinetically irreversible and reversible chemical reactions and to exclude from consideration chemical reactions which remained incomplete in a limited size engine. We present the features of the process of thermal decomposition of liquid and solid fuels which can be used in high-speed aircraft engines.     

Computing stationary distributions in equilibrium and nonequilibrium systems with forward flux sampling

We present a method for computing stationary distributions for activated processes in equilibrium and nonequilibrium systems using forward flux sampling. In this method, the stationary distributions are obtained directly from the rate constant calculations for the forward and backward reactions; there is no need to perform separate calculations for the stationary distribution and the rate constant. We apply the method to the nonequilibrium rare event problem proposed by Maier and Stein, to nucleation in a 2-dimensional Ising system, and to the flipping of a genetic switch.     

Structure, thermodynamics, and liquid-vapor equilibrium of ethanol from molecular-dynamics simulations using nonadditive interactions

We present a molecular-dynamics simulation study of the bulk and liquid-vapor interfacial properties of ethanol using a polarizable force field based on the fluctuating charge (FQ) formalism, as well as the nonpolarizable CHARMM22 force field. Both models are competitive with respect to the prediction of ambient liquid properties such as liquid density, enthalpy of vaporization, dielectric constant, and self-diffusion constants. The polarizable model predicts an average condensed-phase dipole moment of 2.2 D associated with an induced liquid-phase dipole moment of 0.6 D; though qualitatively in agreement with earlier nonadditive models as well as recent Car-Parinello calculations, the current FQ model underestimates the condensed-phase dipole moment. In terms of liquid structure, both models are in agreement with recent neutron-diffraction results of liquid ethanol structure, although the polarizable model predicts the hydroxyl-hydrogen-hydroxyl-hydrogen structure factor in closer agreement with the experimental data. In terms of interfacial properties, both models predict ambient surface tension to within 4% of the experimental value of 22.8 dyncm, while overestimating the surface excess entropy by almost a factor of 2. Both models display the characteristic preferential orientation of interfacial molecules. The polarizable model allows for a monotonic variation of the average molecular dipole moment from the bulk value to that of the vapor phase. Consequently, there is a dramatic difference in the surface potential predicted by the polarizable and nonpolarizable models. The polarizable model estimates a surface potential of -209+/-3 mV, while the nonpolarizable model yields a value of -944+/-10 mV. Finally, based on the vapor-liquid equilibrium simulation data from several temperatures, we estimate the critical properties of both models. As observed with other FQ models for associating fluids (such as water and methanol), and counter to what one would anticipate by modeling more physically the electrostatic response to local environment, the current FQ model underestimates the critical temperature and overestimates the critical density of ethanol; moreover, the FQ model is, in this respect, equivalent to the underlying fixed-charge model. These results further suggest the need to revisit polarizable models in terms of quantitative vapor-liquid equilibrium prediction.     

Experimental investigation of evaporation and condensation in the contact line region of a thin liquid film experiencing small thermal perturbations

Image-analyzing interferometry technique is successfully used to study microscale transport processes related to a curved microfilm on a solid substrate. Digital image processing is used to analyze the images of interference fringes, leading to the evaluation of liquid (heptane) film thickness and curvature profiles at different inclinations on a high refractive index glass surface. The curvature profiles obtained at different inclinations clearly demonstrate that there is a maximum in curvature near the junction of the adsorbed film (of uniform thickness) and the curved film, and then it becomes constant in the thicker portions of the film. The adsorbed film thickness is measured for horizontal as well as inclined positions. Experimentally obtained values of the dispersion constants are compared to those predicted from the Dzyaloshinskii-Lifshitz-Pitaevskii (DLP) theory, and reasonable agreements were obtained. A parameter alpha is defined and experimentally evaluated to quantify the closeness of the system to equilibrium. The nonequilibrium behavior of this parameter alpha is also observed with certain heat input at a particular inclination. A small thermal perturbation is used to force the liquid meniscus to undergo a cycle of alternating condensation and evaporation. High-speed video-microscopy and subsequent image analysis are used for data analysis. The numerical solution of a model that takes into account the balance between the suction and the capillary force is compared with the data to elicit new insights into the evaporation/condensation phenomena and to estimate the interfacial temperature differences for near-equilibrium situations.     

Rheological properties and interfacial slip of a multilayer structure under dynamic shear

Interfacial slip between polymer melt under steady shear has been studied using a simplified multilayer structure. In this investigation, interfacial slip under dynamic shear was studied by calculating the angular displacements of the multilayer structure and its component layers. On the basis of the angular displacements, a slip index was defined to quantify the degree of interfacial slip. A relationship governing the rheological behavior of the multilayer structure under slip and nonslip condition was established. These results were correlated with equations derived from consideration of energy equilibrium in the multilayer structure. Polymer multilayer structures of high‐density polyethylene/polystyrene and liquid crystal polymer(LCP)/poly ethylene naphthalate(PEN) were investigated. Of all the polymers investigated, large interfacial slip was found at LCP/PEN interface under dynamic shear. The high rigidity and alignment along the interface of LCP molecules was believed to prevent chain entanglement in the interfacial layer and therefore promote interfacial slip at the interface. © 2005 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 43: 2683–2693, 2005     

Molecular dynamics simulation of the equilibrium liquid–vapor interphase with solidification

The equilibrium structure of the finite, interphase interfacial region that exists between a liquid film and a bulk vapor is resolved by molecular dynamics simulation. Argon systems are considered for a temperature range that extends below the melting point. Physically consistent procedures are developed to define the boundaries between the interphase and the liquid and vapor phases. The procedures involve counting of neighboring molecules and comparing the results with boundary criteria that permit the boundaries to be precisely established. Two-dimensional radial distribution functions at the liquid and vapor boundaries and within the interphase region demonstrate the physical consistency of the boundary criteria and the state of transition within the region. The method developed for interphase boundary definitions can be extended to nonequilibrium systems. Spatial profiles of macroscopic properties across the interphase region are presented. A number of interfacial thermodynamic properties and profile curve-fit parameters are tabulated, including evaporation/condensation coefficients determined from molecular flux statistics. The evaporation/condensation coefficients away from the melting point compare more favorably with transition state theory than those of previous simulations. Near the melting point, transition theory approximations are less valid and the present results differ from the theory. The effects of film substrate wetting on evaporation/condensation coefficients are also presented.