Molecular dynamics simulation of the diffusion of nanoconfined fluids

A new molecular dynamics simulation technique for simulating fluids in confinement [H. Eslami, F. Mozaffari, J. Moghadasi, F. Müller-Plathe, J. Chem. Phys. 129 (2008) 194702] is employed to simulate the diffusion coefficient of nanoconfined Lennard-Jones fluid. The diffusing fluid is liquid Ar and the confining surfaces are solid Ar fcc (100) surfaces, which are kept frozen during the simulation. In this simulation just the fluid in confinement is simulated at a constant temperature and a constant parallel component of pressure, which is assumed to be equal to the bulk pressure. It is shown that the calculated parallel (to the surfaces) component of the diffusion coefficients depends on the distance between the surfaces (pore size) and shows oscillatory behavior with respect to the intersurface separations. Our results show that on formation of well-organized layers between the surfaces, the parallel diffusion coefficients decrease considerably with respect to the bulk fluid. The effect of pressure on the parallel diffusion coefficients has also been studied. Better organized layers, and hence, lower diffusion coefficients are observed with increasing the pressure.     

Hybrid method coupling fluctuating hydrodynamics and molecular dynamics for the simulation of macromolecules

We present a hybrid computational method for simulating the dynamics of macromolecules in solution which couples a mesoscale solver for the fluctuating hydrodynamics (FH) equations with molecular dynamics to describe the macromolecule. The two models interact through a dissipative Stokesian term first introduced by Ahlrichs and Dunweg [J. Chem. Phys. 111, 8225 (1999)]. We show that our method correctly captures the static and dynamical properties of polymer chains as predicted by the Zimm model. In particular, we show that the static conformations are best described when the ratio sigma/b=0.6, where sigma is the Lennard-Jones length parameter and b is the monomer bond length. We also find that the decay of the Rouse modes' autocorrelation function is better described with an analytical correction suggested by Ahlrichs and Dunweg. Our FH solver permits us to treat the fluid equation of state and transport parameters as direct simulation parameters. The expected independence of the chain dynamics on various choices of fluid equation of state and bulk viscosity is recovered, while excellent agreement is found for the temperature and shear viscosity dependence of center of mass diffusion between simulation results and predictions of the Zimm model. We find that Zimm model approximations start to fail when the Schmidt number Sc < or approximately 30. Finally, we investigate the importance of fluid fluctuations and show that using the preaveraged approximation for the hydrodynamic tensor leads to around 3% error in the diffusion coefficient for a polymer chain when the fluid discretization size is greater than 50 A.     

Molecular dynamics simulation in the grand canonical ensemble

An extended system Hamiltonian is proposed to perform molecular dynamics (MD) simulation in the grand canonical ensemble. The Hamiltonian is similar to the one proposed by Lynch and Pettitt (Lynch and Pettitt, J Chem Phys 1997, 107, 8594), which consists of the kinetic and potential energies for real and fractional particles as well as the kinetic and potential energy terms for material and heat reservoirs interacting with the system. We perform a nonlinear scaling of the potential energy parameters of the fractional particle, as well as its mass to vary the number of particles dynamically. On the basis of the equations of motion derived from this Hamiltonian, an algorithm has been proposed for MD simulation at constant chemical potential. The algorithm has been tested for the ideal gas, for the Lennard-Jones fluid over a wide range of temperatures and densities, and for water. The results for the low-density Lennard-Jones fluid are compared with the predictions from a truncated virial equation of state. In the case of the dense Lennard-Jones fluid and water our predicted results are compared with the results reported using other available methods for the calculation of the chemical potential. The method is also applied to the case of vapor-liquid coexistence point predictions.     

Grand canonical-like molecular dynamics simulations: application to anisotropic mass diffusion in a nanoporous medium

In this work, we describe two grand canonical-like molecular dynamics approaches to investigate mass diffusion phenomenon of a simple Lennard-Jones fluid confined between solid surfaces and in direct contact with reservoirs. In the first method, the density is used as the control variable in the reservoir whereas it is the pressure in the second method. Both methods provide consistent results, however, the constant density approach is the most efficient with respect to the computational time and implementation. Then, employing the constant density approach, we have studied the transient behavior of the diffusion process associated with the migration of one fluid into another one confined between parallel solid walls. Results have shown that the evolution of molar fraction of the invading fluid follows roughly a 1D diffusion model when the solid phase is weakly or moderately adsorbent with a characteristic time increasing when the pore width decreases. However, when the adsorption is high and the pore width small (i.e., below ten molecular sizes), the apparent mass diffusion in the adsorbed layer is reduced compared to that in the center of the slit pore. Hence, this mass diffusion process becomes a two-dimension phenomenon that must take into account an effective mass diffusion coefficient varying locally.     

Molecular theory of thermal conductivity of the Lennard-Jones fluid

In this paper the thermal conductivity of the Lennard-Jones fluid is calculated by applying the combination of the density-fluctuation theory, the modified free volume theory of diffusion, and the generic van der Waals equation of state. A Monte Carlo simulation method is used to compute the equilibrium pair-correlation function necessary for computing the mean free volume and the coefficient in the potential-energy and virial contributions to the thermal conductivity. The theoretical results are compared with our own molecular dynamics simulation results and with those reported in the literature. They agree in good accuracy over wide ranges of density and temperature examined in molecular dynamics simulations. Thus the combined theory represents a molecular theory of thermal conductivity of the Lennard-Jones fluid and by extension simple fluids, which enables us to compute the nonequilibrium quantity by means of the Monte Carlo simulations for the equilibrium pair-correlation function.     

On the kinetics of spinodal decomposition: a numerical solution of a dissipative wave equation

A numerical solution of a recently derived dissipative wave equation governing the kinetics of Spinodal decomposition of a Lennard-Jones fluid is presented. In addition, the results are compared with those of Cahn'S and Abraham's generalized diffusion theories for the case of the early stages of the coarsening process.     

Thermodiffusion in model nanofluids by molecular dynamics simulations

In this work, a new algorithm is proposed to compute single particle (infinite dilution) thermodiffusion using nonequilibrium molecular dynamics simulations through the estimation of the thermophoretic force that applies on a solute particle. This scheme is shown to provide consistent results for model nanofluids in the liquid state (spherical nonmetallic nanoparticles+Lennard-Jones fluid) where it appears that thermodiffusion amplitude, as well as thermal conductivity, decreases with nanoparticle concentration. Then, by changing the nature of the nanoparticle (size, mass, and internal stiffness) and that of the solvent (quality and viscosity), various trends are exhibited. In all cases, the single particle thermodiffusion is positive, i.e., the nanoparticle tends to migrate toward the cold area. The single particle thermal diffusion coefficient is shown to be independent of the size of the nanoparticle (diameter of 0.8-4 nm), whereas it increases with the quality of the solvent and is inversely proportional to the viscosity of the fluid. In addition, this coefficient is shown to be independent of the mass of the nanoparticle and to increase with the stiffness of the nanoparticle internal bonds. Besides, for these configurations, the mass diffusion coefficient behavior appears to be consistent with a Stokes-Einstein-like law.     

On the kinetics of the capillary imbibition of a simple fluid through a designed nanochannel using the molecular dynamics simulation approach

A molecular dynamics (MD) approach was employed to simulate the imbibition of a designed nanopore by a simple fluid (i.e., a Lennard-Jones (LJ) fluid). The length of imbibition as a function of time for various interactions between the LJ fluid and the pore wall was recorded for this system (i.e., the LJ fluid and the nanopore). By and large, the kinetics of imbibition was successfully described by the Lucas-Washburn (LW) equation, although deviation from it was observed in some cases. This lack of agreement is due to the neglect of the dynamic contact angle (DCA) in the LW equation. Two commonly used models (i.e., hydrodynamic and molecular-kinetic (MK) models) were thus employed to calculate the DCA. It is demonstrated that the MK model is able to justify the simulation results in which are not in good agreement with the simple LW equation. However, the hydrodynamic model is not capable of doing that. Further investigation of the MD simulation data revealed an interesting fact that there is a direct relationship between the wall-fluid interaction and the speed of the capillary imbibition. More evidence to support this claim is presented.     

Lattice-Boltzmann simulations of the dynamics of polymer solutions in periodic and confined geometries

A numerical method to simulate the dynamics of polymer solutions in confined geometries has been implemented and tested. The method combines a fluctuating lattice-Boltzmann model of the solvent [Ladd, Phys. Rev. Lett. 70, 1339 (1993)] with a point-particle model of the polymer chains. A friction term couples the monomers to the fluid [Ahlrichs and Dunweg, J. Chem. Phys. 111, 8225 (1999)], providing both the hydrodynamic interactions between the monomers and the correlated random forces. The coupled equations for particles and fluid are solved on an inertial time scale, which proves to be surprisingly simple and efficient, avoiding the costly linear algebra associated with Brownian dynamics. Complex confined geometries can be represented by a straightforward mapping of the boundary surfaces onto a regular three-dimensional grid. The hydrodynamic interactions between monomers are shown to compare well with solutions of the Stokes equations down to distances of the order of the grid spacing. Numerical results are presented for the radius of gyration, end-to-end distance, and diffusion coefficient of an isolated polymer chain, ranging from 16 to 1024 monomers in length. The simulations are in excellent agreement with renormalization group calculations for an excluded volume chain. We show that hydrodynamic interactions in large polymers can be systematically coarse-grained to substantially reduce the computational cost of the simulation. Finally, we examine the effects of confinement and flow on the polymer distribution and diffusion constant in a narrow channel. Our results support the qualitative conclusions of recent Brownian dynamics simulations of confined polymers [Jendrejack et al., J. Chem. Phys. 119, 1165 (2003) and Jendrejack et al., J. Chem. Phys. 120, 2513 (2004)].     

NVU dynamics. II. Comparing to four other dynamics

In the companion paper [T. S. Ingebrigtsen, S. Toxvaerd, O. J. Heilmann, T. B. Schr?der, and J. C. Dyre, "NVU dynamics. I. Geodesic motion on the constant-potential-energy hypersurface," J. Chem. Phys. (in press)] an algorithm was developed for tracing out a geodesic curve on the constant-potential-energy hypersurface. Here, simulations of NVU dynamics are compared to results for four other dynamics, both deterministic and stochastic. First, NVU dynamics is compared to the standard energy-conserving Newtonian NVE dynamics by simulations of the Kob-Andersen binary Lennard-Jones liquid, its WCA version (i.e., with cut-off's at the pair potential minima), and the Lennard-Jones Gaussian liquid. We find identical results for all quantities probed: radial distribution functions, incoherent intermediate scattering functions, and mean-square displacement as function of time. Arguments are presented for the equivalence of NVU and NVE dynamics in the thermodynamic limit; in particular, to leading order in 1∕N these two dynamics give identical time-autocorrelation functions. In the final part of the paper, NVU dynamics is compared to Monte Carlo dynamics, to a diffusive dynamics of small-step random walks on the constant-potential-energy hypersurface, and to Nose?-Hoover NVT dynamics. If time is scaled for the two stochastic dynamics to make single-particle diffusion constants identical to that of NVE dynamics, the simulations show that all five dynamics are equivalent at low temperatures except at short times.     

Mutual diffusion coefficients of heptane isomers in nitrogen: a molecular dynamics study

The accurate knowledge of transport properties of pure and mixture fluids is essential for the design of various chemical and mechanical systems that include fluxes of mass, momentum, and energy. In this study we determine the mutual diffusion coefficients of mixtures composed of heptane isomers and nitrogen using molecular dynamics (MD) simulations with fully atomistic intermolecular potential parameters, in conjunction with the Green-Kubo formula. The computed results were compared with the values obtained using the Chapman-Enskog (C-E) equation with Lennard-Jones (LJ) potential parameters derived from the correlations of state values: MD simulations predict a maximum difference of 6% among isomers while the C-E equation presents that of 3% in the mutual diffusion coefficients in the temperature range 500-1000 K. The comparison of two approaches implies that the corresponding state principle can be applied to the models, which are only weakly affected by the anisotropy of the interaction potentials and the large uncertainty will be included in its application for complex polyatomic molecules. The MD simulations successfully address the pure effects of molecular structure among isomers on mutual diffusion coefficients by revealing that the differences of the total mutual diffusion coefficients for the six mixtures are caused mainly by heptane isomers. The cross interaction potential parameters, collision diameter σ(12), and potential energy well depth ?(12) of heptane isomers and nitrogen mixtures were also computed from the mutual diffusion coefficients.     

Improvement of methods for molecular dynamics integration

An application of symplectic implicit Runge–Kutta (RK ) integration schemes, the s-stage Gauss–Legendre Runge–Kutta (GLRK ) methods of order 2s, for the numerical solution of molecular dynamics (MD ) equation is described. The two-stage fourth-order GLRK method, the implicit midpoint rule, and the three-stage diagonally implicit RK method of order four are studied. The fixed-point iteraction was used for solving the resulting nonlinear system of equations. The algorithms were applied to a complex system of N particles interacting through a Lennard-Jones potential. The proposed symplectic methods for MD integration permit a wide range of time steps, are highly accurate and stable, and are thus suitable for the MD integration. © 1994 John Wiley & Sons, Inc.     

Application of Fokker-Planck-Kramers equation treatment for short-time dynamics of diffusion-controlled reaction in supercritical Lennard-Jones fluids over a wide density range

The validity of a Fokker-Planck-Kramers equation (FPKE) treatment of the rate of diffusion-controlled reaction at short times [K. Ibuki and M. Ueno, J. Chem. Phys. 119, 7054 (2003)] is tested in a supercritical Lennard-Jones fluid over a wide density range by comparing it with the Langevin dynamics and molecular dynamics simulations and other theories. The density n range studied is 0.323n(c)< or =n< or =2.58n(c) and the temperature 1.52T(c), where n(c) and T(c) are the critical density and temperature, respectively. For the rate of bimolecular reactions, the transition between the collision-limited and diffusion-limited regimes is expected to take place in this density range. The simulations show that the rate constant decays with time extensively at high densities, and that the magnitude of decay decreases gradually with decreasing density. The decay profiles of the rate constants obtained by the simulations are reproduced reasonably well by the FPKE treatment in the whole density range studied if a continuous velocity distribution is used in solving the FPKE approximately. If a discontinuous velocity distribution is used instead of the continuous one, the FPKE treatment leads to a rate constant much larger than the simulation results at medium and low densities. The rate constants calculated from the Smoluchowski-Collins-Kimball (SCK) theory based on the diffusion equation are somewhat smaller than the simulation results in medium and low densities when the intrinsic rate constant is chosen to adjust the steady state rate constant in the low density limit to that derived by the kinetic collision theory. The discrepancy is relatively small, so that the SCK theory provides a useful guideline for a qualitative discussion of the density effect on the rate constant.     

Theoretical description of dissipative vibrational dynamics using the density matrix in the state representation

Ultrafast dissipative dynamics of vibrational degrees of freedom in molecular systems in the condensed phase are studied here. Assuming that the total system is separable into a relevant part and a reservoir, the dynamics of the relevant part can be described by means of a reduced statistical density operator. For a weak or intermediate coupling between the relevant part and the reservoir, it is possible to derive a second-order master equation for this operator. Using a representation of the reduced statistical operator in an appropriate molecular basis set, vibrational dynamics in a variety of potential energy surfaces can be studied. In the numerical calculations, we focus on the dissipative dynamics under the influence of external laser fields. In the first example, vibrational wave-packet dynamics and time-resolved pump-probe spectroscopy of molecular systems with nonadiabatically coupled excited-state potential energy surfaces is presented. In the second part, we show how an intense laser field modifies the wave-packet motion onto two radiatively coupled potential energy surfaces. Finally, the controlled preparation of definite vibrational states in a triatomic molecule with infrared laser pulses is considered taking relaxation and dephasing processes into account. © 1996 John Wiley & Sons, Inc.     

Finite system size effects in the interfacial dynamics of binary liquid films

We study the relaxation dynamics of capillary waves in the interface between two confined liquid layers by means of molecular dynamics simulations. We measure the autocorrelations of the interfacial Fourier modes and find that the finite thickness of the liquid layers leads to a marked increase of the relaxation times as compared to the case of fluid layers of infinite depth. The simulation results are in good agreement with a theoretical first-order perturbation derivation, which starts from the overdamped Stokes' equation. The theory also takes into account an interfacial friction, but the difference with no-slip interfacial conditions is small. When the walls are sheared, it is found that the relaxation times of modes perpendicular to the flow are unaffected. Modes along the flow direction are relatively unaffected as long as the equilibrium relaxation time is sufficiently short compared to the rate of deformation. We discuss the consequences for experiments on thin layers and on ultralow surface tension fluids, as well as computer simulations.     

Transient molecular dynamics simulations of viscosity for simple fluids

A transient molecular dynamics (TMD) method has been developed for simulation of fluid viscosity. In this method a sinusoidal velocity profile is instantaneously overlaid onto equilibrated molecular velocities, and the subsequent decay of that velocity profile is observed. The viscosity is obtained by matching in a least-squares sense the analytical solution of the corresponding momentum transport boundary-value problem to the simulated decay of the initial velocity profile. The method was benchmarked by comparing results obtained from the TMD method for a Lennard-Jones fluid with those previously obtained using equilibrium molecular dynamics (EMD) simulations. Two different constitutive models were used in the macroscopic equations to relate the shear rate to the stress. Results using a Newtonian fluid model agree with EMD results at moderate densities but exhibit an increasingly positive error with increasing density at high densities. With the initial velocity profiles used in this study, simulated transient velocities displayed clear viscoelastic behavior at dimensionless densities above 0.7. However, the use of a linear viscoelastic model reproduces the simulated transient velocity behavior well and removes the high-density bias observed in the results obtained under the assumption of Newtonian behavior. The viscosity values obtained using the viscoelastic model are in excellent agreement with the EMD results over virtually the entire fluid domain. For simplicity, the Newtonian fluid model can be used at lower densities and the viscoelastic model at higher densities; the two models give equivalent results at intermediate densities.     

Simulation of reorientation dynamics in bipolar nematic droplets

A two-dimensional model composed of a synthesis of the Leslie-Ericksen continuumtheory of nematics and the Euler-Lagrange equation for surface director motion is used to study the magnetic-induced director reorientation dynamics confined in spherical bipolar droplets with viscoelastic surfaces. The magnetic field is restricted to the droplet axis of symmetry direction. The numerical results indicate that the surface viscosity and anchoring strength must be taken into account to describe accurately director reorientation dynamics in droplets. In addition, the numerical results replicate frequently reported experimental observations on the performance of polymer dispersed liquid crystal films. These observations include the familiar exponential increase followed by saturation in light transmittance as the external applied field increases, and the exponential increase (decrease) followed by saturation as time increases in the on (off) state. Furthermore, this model is able to predict precisely the relationships between the rise and decay times and the external applied field strength, and the fact that the switching field strength is inversely proportional to droplet size.     

Phase transitions of methane using molecular dynamics simulations

Using a short ranged Lennard-Jones interaction and a long ranged electrostatic potential, CH4 under high pressure was modeled. Molecular dynamics simulations on small clusters (108 and 256 molecules) were used to explore the phase diagram. Regarding phase transitions at different temperatures, our numerical findings are consistent with experimental results to a great degree. In addition, the hysteresis effect is displayed in our results.     

Effective stochastic dynamics on a protein folding energy landscape

We present an approach to protein folding kinetics using stochastic reaction-coordinate dynamics, in which the effective drift velocities and diffusion coefficients are determined from microscopic simulation data. The resultant Langevin equation can then be used to directly simulate the folding process. Here, we test this approach by applying it to a toy two-state dynamical system and to a funnellike structure-based (Go-type) model. The folding time predictions agree very well with full simulation results. Therefore, we have in hand a fast numerical tool for calculating the folding kinetic properties, even when full simulations are not feasible. In addition, the local drift and diffusion coefficients provide an alternative way to compute the free energy profile in cases where only local sampling can be achieved.     

Nonequilibrium molecular dynamics simulation to describe the rotation of rigid, low aspect ratio carbon nanotubes in simple shear flow

In this paper, the rotation of short carbon nanotubes in simple shear liquid argon flow was investigated by nonequilibrium molecular dynamics (MD) simulation. In their simulations, nanotubes were described as rigid cylinders of carbon atoms. Lennard-Jones potential was employed to represent both argon-argon and argon-carbon interactions. Results show that time period of a nanotube as calculated from MD simulations is longer than what would be calculated from Jeffery's equation based on the aspect ratio of the cylinder. The difference is much higher at low shear rates and for small aspect ratios. Results also reveal that adding caps to an open-ended nanotube speeds up its rotation.