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.     

Viscosity of heptane-toluene mixtures. Comparison of molecular dynamics and group contribution methods

Three methods of molecular dynamics simulation [Green–Kubo (G–K), non-equilibrium molecular dynamics (NEMD) and reversed non-equilibrium molecular dynamics (RNEMD)], and two group contribution methods [UNIFAC–VISCO and Grunberg–Nissan (G–N)] were used to calculate the viscosity of mixtures of n-heptane and toluene (known as heptol). The results obtained for the viscosity and density of heptol were compared with reported experimental data, and the advantages and disadvantages of each method are discussed. Overall, the five methods showed good agreement between calculated and experimental viscosities. In all cases, the deviation was lower than 9%. It was found that, as the concentration of toluene increases, the deviation of the density of the mixture (as calculated with molecular dynamics methods) also increases, which directly affects the viscosity result obtained. Among the molecular simulation techniques evaluated here, G–K produced the best results, and represents the optimal balance between quality of result and time required for simulation. The NEMD method produced acceptable results for the viscosity of the system but required more simulation time as well as the determination of an appropriate shear rate. The RNEMD method was fast and eliminated the need to determine a set of values for shear rate, but introduced large fluctuations in measurements of shear rate and viscosity. The two group contribution methods were accurate and fast when used to calculate viscosity, but require knowledge of the viscosity of the pure compounds, which is a serious limitation for applications in complex multicomponent systems.     

The Rheology of 6-Propyl Duodecane and 5-Dibutyl Nonane by Non-Equilibrium Molecular Dynamics Simulations

Abstract

In recent papers, we reported non-equilibrium molecular dynamics (NEMD) simulations of planar Couette flow for liquid n- and i-butane, and liquid n-decane and 4-propyl heptane, using two collapsed atom models and an atomistically detailed model. It was found that the collapsed atom models predict the viscosities of the n-butane and n-decane quite well, and that the atomistically detailed model does not yield quantitative agreement with the viscosity of the n-alkanes or the branched alkanes, but it does have the one positive feature that the calculated viscosities of the branched alkanes are higher than these of the n-alkanes. In the present paper, we report results of NEMD simulations of planar Couette flow for liquid 6-propyl duodecane and 5-dibutyl nonane at 296 K and 0.782 g/cc, using an expanded collapsed atom model for simplicity. The strain rate dependent viscosity shows three different regions—two shear thinning ones and a Newtonian one. The slopes of the log-log plot for the branched-chain alkanes at high strain rate are different from those at intermediate strain rate, which is characterized as a rheological behavior of branched-chain alkanes. The Newtonian viscosity of the branched-chain alkanes can be extrapolated from the plateau value of the strain rate dependent viscosity at low strain rate as for straight-chain alkanes [J. Chem. Phys., 105, 1214 (1996)]. The results indicate that more-branched alkanes have a larger viscosity than less-branched C₁₇ alkanes.     

Accurate determination of normal stress differences via transient-time correlation function – non-equilibrium molecular dynamics (TTCF–NEMD) simulations

When a non-Newtonian liquid undergoes shear flow, differences between the values taken by the diagonal elements of the pressure tensor arise. These differences, known as normal stress differences, are notoriously difficult to compute accurately. In this work, we propose to apply the transient-time correlation function (TTCF) formalism together with non-equilibrium molecular dynamics (NEMD) simulations to determine normal stress differences, as well as the corresponding normal stress coefficients, in an atomic fluid undergoing shear flow. Our results shed light on the dependence of normal stress differences and coefficients on the applied shear.     

Long time dynamics of Met-enkephalin: comparison of explicit and implicit solvent models

Met-enkephalin is one of the smallest opiate peptides. Yet, its dynamical structure and receptor docking mechanism are still not well understood. The conformational dynamics of this neuron peptide in liquid water are studied here by using all-atom molecular dynamics (MD) and implicit water Langevin dynamics (LD) simulations with AMBER potential functions and the three-site transferable intermolecular potential (TIP3P) model for water. To achieve the same simulation length in physical time, the full MD simulations require 200 times as much CPU time as the implicit water LD simulations. The solvent hydrophobicity and dielectric behavior are treated in the implicit solvent LD simulations by using a macroscopic solvation potential, a single dielectric constant, and atomic friction coefficients computed using the accessible surface area method with the TIP3P model water viscosity as determined here from MD simulations for pure TIP3P water. Both the local and the global dynamics obtained from the implicit solvent LD simulations agree very well with those from the explicit solvent MD simulations. The simulations provide insights into the conformational restrictions that are associated with the bioactivity of the opiate peptide dermorphin for the delta-receptor.     

Coarse-grained modelling of urea-adamantyl functionalised poly(propylene imine) dendrimers

To investigate the behaviour of poly(propylene imine) dendrimers – and urea–adamantyl functionalised ones – in solution using molecular dynamics simulations, we developed a coarse-grained model to tackle the relatively large system sizes and time scales needed. Harmonic bond and angle potentials were derived from atomistic simulations using an iterative Boltzmann inversion scheme, modified to incorporate Gaussian fits of the bond and angle distributions. With the coarse-grained model and accompanying force field simulations of generations 1–7 of both dendrimer types in water were performed. They compare favourably with atomistic simulations and experimental results on the basis of size, shape, monomer density, spacer back-folding and atomic form factor measurements. These results show that the structural dynamics of these dendrimers originate from flexible chains constrained by configurational and spatial requirements. Large dendrimers are more rigid and spherical, while small ones are flexible, alternatively rod-like and globular.     

Calculating the bulk modulus for a lipid bilayer with nonequilibrium molecular dynamics simulation

Nonequilibrium molecular dynamics (NEMD) computer simulations are used to calculated the bulk modulus for a dimyristoylphosphatidylcholine bilayer. A methodology is developed whereby NEMD can be effectively used to calculate material properties for complex systems that undergo long time-scale conformational changes. It is found that the bulk modulus upon expansion from a zero stress state agrees well with experimental estimates. However, it is also found that the modulus upon contraction from a zero stress state is larger. From a molecular perspective, it is possible to explain this phenomena by examining the molecular origins of the pressure response. The finding that the two moduli are not equal upon compression and expansion is in apparent contradiction to osmotic stress experiments where the area modulus was found to be the same upon expansion and contraction. This issue is addressed.     

Molecular dynamics studies in nanoscale liquid structures: geometry and thermal effects on nanojet development

Conventional macroscopic jet theory relies heavily on experimental correlations which cannot be easily extended to the nanoscale regime. Moreover, the fluid dynamic effects at small length scales and their contribution to the development of nanoscale liquid structures are fundamentally different from their macroscopic counterparts. This coupled with the high spatial and temporal resolution requirements at nanoscale domains make molecular dynamics (MD) an excellent tool for studying such structures. In this study, the formation and breakup of nanojets (NJs) developing from high pressure into vacuum is investigated using MD based on non-Hamiltonian formulations. By ejecting the equilibrated argon atoms through various nozzle geometries and diameters, nanoscale jet flows were generated. The dependence of the jet structure on nozzle geometry and diameter is studied. The influence of geometry on NJ formation is also studied along with issues involved in the equilibration and thermostat coupling parameter. Various thermostats are compared to understand the role they play in MD simulations of liquid nanostructures. Tuning of the thermostat coupling parameter has also been discussed. The jet breakup phenomenon is analysed and a comparative study, vis-à-vis, well-established continuum and stochastic models, is attempted.     

Computer Simulation of Liquid Methanol II. System Size Effects

Abstract

A series of molecular dynamics simulations of liquid methanol has been carried out on a supernode transputer array. Four system sizes from 125 to 512 molecules have been considered, in order to study the effect of system size on the calculated structural, orientational and dynamic properties. The dielectric constant and the dielectric relaxation time are compared with experimental data.     

Intermonolayer friction and surface shear viscosity of lipid bilayer membranes

The flow behavior of lipid bilayer membranes is characterized by a surface viscosity for in-plane shear deformations, and an intermonolayer friction coefficient for slip between the two leaflets of the bilayer. Both properties have been studied for a variety of coarse-grained double-tailed model lipids, using equilibrium and nonequilibrium molecular dynamics simulations. For lipids with two identical tails, the surface shear viscosity rises rapidly with tail length, while the intermonolayer friction coefficient is less sensitive to the tail length. Interdigitation of lipid tails across the bilayer midsurface, as observed for lipids with two distinct tails, strongly enhances the intermonolayer friction coefficient, but hardly affects the surface shear viscosity. The simulation results are compared against the available experimental data.     

The dynamic stress responses to area change in planar lipid bilayer membranes

The viscoelastic properties of planar phospholipid (dimyristoylphosphatidylcholine) bilayer membranes at 308 K are studied, many of them for the first time, using the nonequilibrium molecular dynamics simulation (NEMD) method for membrane area change. First, we present a unified formulation of the intrinsic three-dimensional (3D) and apparent in-plane viscoelastic moduli associated with area change based on the constitutive relations for a uniaxial system. The NEMD simulations of oscillatory area change process are then used to obtain the frequency-domain moduli. In the 4-250 GHz range, the intrinsic 3D elastic moduli of 20-27 kbar and viscous moduli of 0.2-9 kbar are found with anisotropy and monotonic frequency dispersion. In contrast, the apparent in-plane elastic moduli (1-9 kbar) are much smaller than, and the viscous moduli (2-6 kbar) comparable to, their 3D counterparts, due to the interplay between the lateral and normal relaxations. The time-domain relaxation functions, separately obtained by applying stepwise strains, can be fit by 4-6 exponential decay modes spanning subpicosecond to nanosecond timescale and are consistent with the frequency-domain results. From NEMD with varying strain amplitude, the linear constitutive model is shown to be valid up to 6 and 20% area change for the intrinsic 3D elastic and viscous responses, respectively, and up to 20% area change for the apparent in-plane viscoelasticity. Inclusion of a gramicidin A dimer (approximately 1 mol %) yields similar response properties with possibly smaller (<10%) viscous moduli. Our results agree well with available data from ultrasonic experiments, and demonstrate that the third dimension (thickness) of the planar lipid bilayer is integral to the in-plane viscoelasticity.     

Shear and extensional rheology of cellulose/ionic liquid solutions

In this study, we characterize the shear and extensional rheology of dilute to semidilute solutions of cellulose in the ionic liquid 1-ethyl-3-methylimidazolium acetate (EMIAc). In steady shear flow, the semidilute solutions exhibit shear thinning, and the high-frequency complex modulus measured in small amplitude oscillatory shear flow exhibits the characteristic scaling expected for solutions of semiflexible chains. Flow curves of the steady shear viscosity plotted against shear rate closely follow the frequency dependence of the complex viscosity acquired using oscillatory shear, thus satisfying the empirical Cox-Merz rule. We use capillary thinning rheometry (CaBER) to characterize the relaxation times and apparent extensional viscosities of the semidilute cellulose solutions in a uniaxial extensional flow that mimics the dynamics encountered in the spin-line during fiber spinning processes. The apparent extensional viscosity and characteristic relaxation times of the semidilute cellulose/EMIAc solutions increase dramatically as the solutions enter the entangled concentration regime at which fiber spinning becomes viable.     

Leukocytes rolling and recruitment by endothelial cells: hemorheological experiments and numerical simulations

The recruitment of leukocytes from the blood stream and their subsequent adhesion to endothelial walls are essential stages to the immune response system during inflammation. The precise dynamic mechanisms by which molecular mediators facilitate leukocyte arrests are still unknown. In this study combined experimental results and computer simulations are used to investigate localized hydrodynamics of individual and collective behavior of clusters of leukocytes. Leukocyte-endothelial cell interactions in post-capillary venules of Wistar rats cremaster muscle were monitored by intravital microscopy. From these experiments the hemorheologic and hemodynamical measured parameters were used in time dependent three-dimensional computer simulations, using a mesoscopic lattice Boltzmann flow solver for shear thinning fluids. The dynamics of leukocyte clusters under generalized Newtonian blood flow with shear thinning viscosity was computed and discussed. In this paper we present quantified distributions of velocity and shear stress on the surface of leukocytes and near vessel wall attachment points. We have observed one region of maximum shear stress and two regions of minimum shear stress on the surface of leukocytes close to the endothelial wall. We verified that the collective hydrodynamic behavior of the cluster of recruited leukocytes establishes a strong motive for additional leukocyte recruitment. It was found that the lattice Boltzmann solver used here is fully adaptive to the measured experimental parameters. This study suggests that the influence of the leukocytes rolling on the increase of the endothelial wall shear stress may support the activation of more signalling mediators during inflammation.     

Molecular dynamics of myoglobin at 298 degrees K. Results from a 300-ps computer simulation.

We have carried out a very long (300 ps) molecular dynamics simulation of the protein myoglobin. This trajectory is approximately three times longer than the longest previous molecular dynamics simulation of a protein, and ten times longer than protein simulations of comparable size (1,423 atoms in our model). Here we report results from this long simulation concerning the average structure, the mean square fluctuations of atoms about the average structure, and the nuclear magnetic resonance order parameters for various groups in myoglobin. The results demonstrate that the average coordinates change very slowly during the simulation. The relative atomic mobilities are well described by the simulation. For both the mean square atomic fluctuations and the order parameters, however, there are significant quantitative differences when values calculated using shorter portions of the trajectory are compared with results obtained for the entire 300-ps simulation. The implications of this result for obtaining converged properties from protein molecular dynamics simulations for comparison with experiment are discussed.     

Typical trajectories of coupled degrade-and-fire oscillators: from dispersed populations to massive clustering

We consider the dynamics of a piecewise affine system of degrade-and-fire oscillators with global repressive interaction, inspired by experiments on synchronization in colonies of bacteria-embedded genetic circuits. Due to global coupling, if any two oscillators happen to be in the same state at some time, they remain in sync at all subsequent times; thus clusters of synchronized oscillators cannot shrink as a result of the dynamics. Assuming that the system is initiated from random initial configurations of fully dispersed populations (no clusters), we estimate asymptotic cluster sizes as a function of the coupling strength. A sharp transition is proved to exist that separates a weak coupling regime of unclustered populations from a strong coupling phase where clusters of extensive size are formed. Each phenomena occurs with full probability in the thermodynamics limit. Moreover, the maximum number of asymptotic clusters is known to diverge linearly in this limit. In contrast, we show that with positive probability, the number of asymptotic clusters remains bounded, provided that the coupling strength is sufficiently large.     

Predicted melt curve and liquid-state transport properties of TATB from molecular dynamics simulations

The melt curve and the liquid-state transport properties shear viscosity, self-diffusion coefficient and thermal conductivity of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) were predicted using all-atom molecular dynamics simulations. The TATB melt curve was obtained using solid–liquid coexistence simulations and is in good accord with the Simon–Glatzel equation. The temperature dependencies of the shear viscosity and self-diffusion coefficient are predicted to obey Arrhenius behaviour for pressures up to P = 20 kbar. The thermal conductivity has a linear temperature dependence for P < 15 kbar and a linear density (ρ) dependence for ρ > 1200 kg m^?3. At similar densities the shear viscosity of liquid TATB is close to the predictions for liquid nitromethane [58] but lower than the predictions for liquid HMX [24] and RDX [59]. The self-diffusion coefficient for TATB is predicted to be higher than predictions for nitromethane, HMX and RDX at similar densities. The conductivity of TATB is ≈20% greater than the conductivity of liquid HMX at a given density.     

Computer Simulation of Isothermal Mass Transport in Graphite Slit Pores

Abstract

Results are presented from a simulation study of the mass transport of oxygen and nitrogen through graphite slit pores. The work is motivated by an attempt to understand the molecular origins of the kinetic selectivity displayed when air is separated into its major components using pressure swing adsorption. A combination of non-equilibrium molecular dynamics (NEMD), equilibrium molecular dynamics (EMD) and grand canonical Monte Carlo methods has been employed in our study to extract the maximum information. Transport diffusivities, self-diffusivities, permeabilities and Darken thermodynamic factors have been calculated as a function of pore width and temperature for pure component oxygen and nitrogen. In addition, new EMD simulation data for an 80:20 mixture of nitrogen and oxygen is reported, including a direct calculation of the Stefan-Maxwell coefficients. The results are discussed in terms of the oxygen selectivity and the possible mechanisms, which increase or decrease this quantity.

We find that the pore width behaviour of the diffusion coefficients consists of three distinct regimes: a regime at larger pore widths in which single component diffusion coefficients are largely independent of pore width, an optimum pore width at which both diffusivities increase substantially but the slit pore is selective towards nitrogen, and a regime at very low pore widths at which the diffusivities decrease sharply, but the slits are selective towards oxygen. The mechanism behind each of these regimes is discussed in terms of “entropic” effects and potential barrier heights.

We have also found that permeability selectivity is substantially reduced in a mixture of the two gases with a composition similar to that of air. Cross diffusion coefficients in the mixture have been calculated and shown to be non-negligible.     

On the relative importance of rheology for image-based CFD models of the carotid bifurcation

BACKGROUND: Patient-specific computational fluid dynamics (CFD) models derived from medical images often require simplifying assumptions to render the simulations conceptually or computationally tractable. In this study, we investigated the sensitivity of image-based CFD models of the carotid bifurcation to assumptions regarding the blood rheology. METHOD OF APPROACH: CFD simulations of three different patient-specific models were carried out assuming: a reference high-shear Newtonian viscosity, two different non-Newtonian (shear-thinning) rheology models, and Newtonian viscosities based on characteristic shear rates or, equivalently, assumed hematocrits. Sensitivity of wall shear stress (WSS) and oscillatory shear index (OSI) were contextualized with respect to the reproducibility of the reconstructed geometry, and to assumptions regarding the inlet boundary conditions. RESULTS: Sensitivity of WSS to the various rheological assumptions was roughly 1.0 dyn/cm(2) or 8%, nearly seven times less than that due to geometric uncertainty (6.7 dyn/cm(2) or 47%), and on the order of that due to inlet boundary condition assumptions. Similar trends were observed regarding OSI sensitivity. Rescaling the Newtonian viscosity based on time-averaged inlet shear rate served to approximate reasonably, if overestimate slightly, non-Newtonian behavior. CONCLUSIONS: For image-based CFD simulations of the normal carotid bifurcation, the assumption of constant viscosity at a nominal hematocrit is reasonable in light of currently available levels of geometric precision, thus serving to obviate the need to acquire patient-specific rheological data.     

Self diffusion of argon in flexible,single wall,carbon nanotubes

The high-throughput Condor environment now allows many simulations to be performed on related systems, whether the focus is on improving the statistics or on broadening the range of conditions under which these simulations run. We illustrate the scope of the approach by using equilibrium molecular dynamics (EMD) to calculate self-diffusivities of argon atoms diffusing through single wall carbon nanotubes (SWNT). The diameters of the tubes and their helicities were varied and different argon loadings were studied. We also considered the effect of the rigidity/flexibility of the tube on the diffusivity. We found that the helicity and flexibility of the tubes have almost no noticeable influences. The size of the pore had a small effect, but the diffusivity depended essentially on the fluid loading.