Invariant dynamics in a united-atom model of an ionic liquid.

We study a united-atom model of the ionic liquid 1-butyl-1-methylpyrrolidinium bis(trifluoromethyl)sulfonylamide to determine to what extent there exist curves in the phase diagram along which the microscopic dynamics are invariant when expressed in dimensionless, or reduced, form. The initial identification of these curves, termed isodynes, is made by noting that contours of reduced shear viscosity and reduced self-diffusion coefficient coincide to a good approximation. Choosing specifically the contours of reduced viscosity as nominal isodynes, further simulations were carried out for state points on these, and other aspects of dynamics were investigated to study their degree of invariance. These include the mean-squared displacement, shear-stress autocorrelation function, and various rotational correlation functions. These were invariant to a good approximation, with the main exception being rotations of the anion about its long axis. The dynamical features that are invariant have in common that they are aspects that would be relevant for a coarse-grained description of the system; specifically, removing the most microscopic degrees of freedom in principle leads to a simplification of the potential energy landscape, which allows for the existence of isodynes.

A comprehensive study of the thermal conductivity of the hard sphere fluid and solid by molecular dynamics simulation.

This work reports a new set of hard sphere (HS) thermal conductivity coefficient, λ, data obtained by Molecular Dynamics (MD) computer simulation, over a density range covering the dilute fluid to near the close-packed solid, and for a large number of particles (up to N = 13 1072) and long simulation times. The N-dependence of the thermal conductivity is shown to be proportional to N-2/3 to a good approximation over a wide range of system sizes, which enabled λ values in the thermodynamic limit to be predicted accurately. The fluid and solid λ can be represented well by the Enskog theory (ET) formula, λE, times a density-dependent correction term, which is close to unity for the fluid and practically constant for the solid. The convergence of the MD λ data back towards ET in the metastable fluid starts just above the freezing density. For the HS solid and dense fluid it was found that the thermal conductivity is nearly linear in pressure, as has been observed experimentally for a number of solids. Simple excess entropy scaling over the higher density fluid phase region was found, and Rosenfeld's exponential relationship can be fitted to the simulation data for the solid to a high degree of accuracy. The simulation analysis has revealed a number of new trends in the behaviour of the HS thermal conductivity which could be useful in building more accurate models for heat conduction in experimental systems.

Thermodynamic and dynamical properties of the hard sphere system revisited by molecular dynamics simulation.

Revised thermodynamic and dynamical properties of the hard sphere (HS) system are obtained from extensive molecular dynamics calculations carried out with large system sizes (number of particles, N) and long times. Accurate formulas for the compressibility factor of the HS solid and fluid branches are proposed, which represent the metastable region and take into account its divergence at close packing. Some basic second-order thermodynamic properties are obtained and a maximum in some of their derivatives in the metastable fluid region is found. The thermodynamic parameters associated with the melting-freezing transition have been determined to four digit accuracy, which generates accurate new values for the coexistence properties of the HS system. For the self-diffusion coefficient, D, it is shown that relatively large systems (N > 104) are required to achieve an accurate linear extrapolation of D to the infinite size limit with a D vs. N-1/3 plot. Moreover, it is found that there is a density dependence of the value of the slope in the linear regime. The density dependent correction becomes practically insignificant at higher densities and the hydrodynamic formula found in the literature is still accurate. However, with decreasing density the density dependence of the size correction cannot be neglected, which indicates that other sources of N-dependence, apart from those derived on purely hydrodynamic grounds, may also be important (and as yet unaccounted for). A detailed analytic representation of the density dependence of the HS self-diffusion coefficient and the HS viscosity, η, is given. It is shown that the HS viscosity near freezing and in the metastable region can be described well by the Krieger-Dougherty equation. Both D and η start to scale at high densities and in the metastable region in such a way that Dηp = const, where p ≃ 0.97, and D → 0 and η → ∞ at a packing fraction of 0.58, which coincides with some previous predictions of the HS glass transition density.

Revisiting the Stokes-Einstein relation without a hydrodynamic diameter.

We present diffusion coefficient and shear viscosity data for the Lennard-Jones fluid along nine isochores above the critical density, each involving a temperature variation of roughly two orders of magnitude. The data are analyzed with respect to the Stokes-Einstein (SE) relation, which breaks down gradually at high temperatures. This is rationalized in terms of the fact that the reduced diffusion coefficient D ̃ and the reduced viscosity η ̃ are both constant along the system's lines of constant excess entropy (the isomorphs). As a consequence, D ̃ η ̃ is a function of T/T Ref(ρ) in which T is the temperature, ρ is the density, and T Ref(ρ) is the temperature as a function of the density along a reference isomorph. This allows one to successfully predict the viscosity from the diffusion coefficient in the studied region of the thermodynamic phase diagram.

Chemical potential of a test hard sphere of variable size in hard-sphere fluid mixtures.

A detailed comparison between the Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) equation of state of hard-sphere mixtures is made with Molecular Dynamics (MD) simulations of the same compositions. The Labík and Smith simulation technique [S. Labík and W. R. Smith, Mol. Simul. 12, 23-31 (1994)] was used to implement the Widom particle insertion method to calculate the excess chemical potential, ßµ0ex, of a test particle of variable diameter, σ0, immersed in a hard-sphere fluid mixture with different compositions and values of the packing fraction, η. Use is made of the fact that the only polynomial representation of ßµ0ex which is consistent with the limits σ0 → 0 and σ0 → ∞ has to be of the cubic form, i.e., c0(η)+c¯1(η)σ0/M1+c¯2(η)(σ0/M1)2+c¯3(η)(σ0/M1)3, where M1 is the first moment of the distribution. The first two coefficients, c0(η) and c¯1(η), are known analytically, while c¯2(η) and c¯3(η) were obtained by fitting the MD data to this expression. This in turn provides a method to determine the excess free energy per particle, ßaex, in terms of c¯2, c¯3, and the compressibility factor, Z. Very good agreement between the BMCSL formulas and the MD data is found for ßµ0ex, Z, and ßaex for binary mixtures and continuous particle size distributions with the top-hat analytic form. However, the BMCSL theory typically slightly underestimates the simulation values, especially for Z, differences which the Boublík-Carnahan-Starling-Kolafa formulas and an interpolation between two Percus-Yevick routes capture well in different ranges of the system parameter space.

Communication: Simple liquids' high-density viscosity.

This paper argues that the viscosity of simple fluids at densities above that of the triple point is a specific function of temperature relative to the freezing temperature at the density in question. The proposed viscosity expression, which is arrived at in part by reference to the isomorph theory of systems with hidden scale invariance, describes computer simulations of the Lennard-Jones system as well as argon and methane experimental data and simulation results for an effective-pair-potential model of liquid sodium.

Chemical potential of a test hard sphere of variable size in a hard-sphere fluid.

The Labík and Smith Monte Carlo simulation technique to implement the Widom particle insertion method is applied using Molecular Dynamics (MD) instead to calculate numerically the insertion probability, P0(η,σ0), of tracer hard-sphere (HS) particles of different diameters, σ0, in a host HS fluid of diameter σ and packing fraction, η, up to 0.5. It is shown analytically that the only polynomial representation of -ln⁡P0(η,σ0) consistent with the limits σ0→0 and σ0→∞ has necessarily a cubic form, c0(η)+c1(η)σ0/σ+c2(η)(σ0/σ)2+c3(η)(σ0/σ)3. Our MD data for -ln⁡P0(η,σ0) are fitted to such a cubic polynomial and the functions c0(η) and c1(η) are found to be statistically indistinguishable from their exact solution forms. Similarly, c2(η) and c3(η) agree very well with the Boublík-Mansoori-Carnahan-Starling-Leland and Boublík-Carnahan-Starling-Kolafa formulas. The cubic polynomial is extrapolated (high density) or interpolated (low density) to obtain the chemical potential of the host fluid, or σ0→σ, as ßµex=c0+c1+c2+c3. Excellent agreement between the Carnahan-Starling and Carnahan-Starling-Kolafa theories with our MD data is evident.

A molecular dynamics study of CaCO3 nanoparticles in a hydrophobic solvent with a stearate co-surfactant.

Stearates containing overbased detergent nanoparticles (NPs) are used as acid neutralising additives in automotive and marine engine oils. Molecular dynamics (MD) simulations of the self-assembly of calcium carbonate, calcium stearate as a co-surfactant and stabilising surfactants of such NPs in a model explicit molecular hydrophobic solvent have been carried out using a methodology described first by Bodnarchuk et al. [J. Phys. Chem. C, 2014, 118, 21092]. The cores and particles as a whole become more elongated with stearate, and the surfactant molecules are more spaced out in this geometry than in their stearate-free counterparts. The rod dimensions are found to be largely independent of the surfactant type for a given amount of CaCO3. The corresponding particles without stearate were more spherical, the precise shape depending to a greater extent on the chemical architecture of the surfactant molecule. The rod-shaped stearate containing nanoparticles penetrated a model water droplet to a greater depth than the corresponding near-spherical particle, which is possibly facilitated by the dissociation of nanoparticle surfactant molecules onto the surface of the water in this process. These simulations are the first to corroborate the nanoparticle-water penetration mechanism proposed previously by experimental groups investigating the NP acid neutralisation characteristics.

Revised Enskog theory and molecular dynamics simulations of the viscosities and thermal conductivity of the hard-sphere fluid and crystal.

Hard-sphere (HS) shear, longitudinal, cross, and bulk viscosities and the thermal conductivity are obtained by molecular dynamics (MD) simulations, covering the entire density range from the dilute fluid to the solid crystal near close-packing. The transport coefficient data for the HS crystal are largely new and display, unlike for the fluid, a surprisingly simple behavior in that they can be represented well by a simple function of the density compressibility factor. In contrast to the other four transport coefficients (which diverge), the bulk viscosity in the solid is quite small and decreases rapidly with increasing density, tending to zero in the close-packed limit. The so-called cross viscosity exhibits a different behavior to the other viscosities, in being negative over the entire solid range, and changes sign from negative to positive on increasing the density in the fluid phase. The extent to which the viscosity tensor and thermal conductivity of the HS crystal can be represented by revised Enskog theory (RET) is investigated. The RET expressions are sums of an instantaneous (I), a kinetic (K), and a so-called α part. The I part of the transport coefficients evaluated directly by MD are statistically indistinguishable from those of the corresponding kinetic theory (Enskog and RET) expressions. For the K part the integral over the spatial two-particle distribution function at contact was determined and the α part was estimated using the direct correlation function and density functional theory approximations. All three parts were determined in this work which allowed the accuracy of RET for solid systems to be assessed rigorously. It is found that in the case of the thermal conductivity the predictions of RET are in excellent agreement with the MD results. Also, for the shear viscosity the agreement over the entire solid phase is quite good but is considerably worse for the three remaining viscosities in the solid phase.

Stationary temperature profiles in a liquid nanochannel: comparisons between molecular-dynamics simulation and classical hydrostatics.

We compare the results of three-dimensional molecular-dynamics (MD) simulations of a Lennard-Jones (LJ) liquid with a hydrostatic (HS) solution of a high temperature liquid channel which is surrounded by a fluid at lower temperature. The maximum temperature gradient, dT/dx , between the two temperature regions ranged from infinity (step function) to dT/dx=0.1 (in the usual LJ units). Because the systems were in stationary-nonequilibrium states with no fluid flow, both MD simulation and the HS solution gave flat profiles for the normal pressure in all temperature-gradient cases. However, the other quantities showed differences between the two methods. The MD-derived density was found to oscillate over the length of ca. 8 LJ particle diameters from the boundary plane in the system with the infinite temperature gradient, while the HS-derived density showed simply a stepwise profile. The MD simulation also showed another anomaly near the boundary in potential energy. We have found systems in which the HS treatment works well and those where the HS approach breaks down, and therefore established the minimum length scale for the HS treatment to be valid. We also compare the kinetic temperature and the configurational temperature in these systems, and show that these can differ in the transition zone between the two temperatures.

First derivative of the hard-sphere radial distribution function at contact.

Molecular dynamics simulations have been carried out of the radial distribution function of the hard sphere fluid for a range of densities in the equilibrium fluid and just into the metastable region. The first derivative of the hard-sphere radial distribution function at contact was computed and its density dependence fitted to a simple analytic form. Comparisons were made with semi-empirical formulae from the literature, and of these the formula proposed by Tao et al (1992 Phys. Rev. A 46 8007) was found to be in best agreement with the simulation data, although it slightly underestimates the derivative at the higher packing fractions in excess of about 0.45. Close to contact, within a few per cent of the particle diameter, the radial distribution function can be represented well by a second order polynomial. An exponential function, which has some useful analytic features, can also be applied in this region.

Comparisons between molecular dynamics and hydrodynamics treatment of nonstationary thermal processes in a liquid.

Molecular dynamics (MD) and Navier-Stokes hydrodynamics have been performed to model thermal relaxation processes arising from an initially established nonequilibrium stationary state. A nanoscale two-layer Lennard-Jones (LJ) liquid system was constructed in which the two parts were initially at a different temperature, with a narrow transitional zone between the two layers that was spatially linear in temperature. The highest-temperature layer had widths of five or 20 LJ particle diameters. The hydrodynamics model used parametrized MD-derived transport coefficients and the LJ equation of state as input functions. The temporal and spatial temperature and density profiles produced by the two methods show good agreement, indicating that a hydrodynamics description is reliable even for nonstationary phenomena down to the scale of a few molecular diameters. We found that at certain locations the Navier-Stokes solution predicted that the pressure and temperature profiles relaxed in a damped oscillatory manner, which we could discern despite the fluctuations in the MD data.

Traction and nonequilibrium phase behavior of confined sheared liquids at high pressure.

Nonequilibrium molecular dynamics simulations of confined model liquids under pressure and sheared by the relative sliding of the boundary walls have been carried out. The relationship between the time-dependent traction coefficient, µ(t), and the state of internal structure of the film is followed from commencement of shear for various control parameters, such as applied load, global shear rate, and solid-liquid atom interaction parameters. Phase diagrams, velocity and temperature profiles, and traction coefficient diagrams are analyzed for pure Lennard-Jones (LJ) liquids and a binary LJ mixture. A single component LJ liquid is found to form semicrystalline arrangements with high-traction coefficients, and stick-slip behavior is observed for high pressures and low-shear velocities, which is shown to involve periodic deformation and stress release of the wall atoms and slip in the solid-liquid boundary region. A binary mixture, which discourages crystallization, gives a more classical tribological response with the larger atoms preferentially adsorbing commensurate with the wall. The results obtained are analyzed in the context of tribology: the binary mixture behaves like a typical lubricant, whereas the monatomic system behaves like a traction fluid. It is discussed how this type of simulation can give insights on the tribological behavior of realistic systems.

Equation of state and structural properties of the Weeks-Chandler-Andersen fluid.

Molecular dynamics simulations have been carried out for the equation of state and percolation properties of the Weeks-Chandler-Andersen (WCA) system in its fluid phase as functions of density and temperature. The compressibility factor Z collapses well for the various isotherms, using an effective particle diameter for the WCA particle which is (in the usual WCA reduced units) sigma(e)=2(16)(1+T)(16), where T is the temperature. A corresponding "effective" packing fraction is zeta(e)=pisigma(e) (3)N6V, for N particles in volume V, which therefore scales out the effects of temperature. Using zeta(e) the simulation derived Z can be fitted to a simple analytic form which is similar to the Carnahan-Starling hard sphere equation of state and which is valid at all temperatures and densities where the WCA fluid is thermodynamically stable. The data, however, are not scalable onto the hard sphere equation of state for the complete packing fraction range. We explored the continuum percolation behavior of the WCA fluids. The percolation distance sigma(p) for the various states collapses well onto a single curve when plotted as sigma(p)sigma(e) against zeta(e). The ratio sigma(p)sigma(e) exhibits a monotonic decrease with increasing zeta(e) between the percolation line for permeable spheres and the glass transition limit, where sigma(p)sigma(e) approximately 1. The percolation packing fraction was calculated as a function of effective packing fraction and fitted to an empirical expression. The local coordination number at the percolation threshold showed a transition between the soft core and hard core limits from ca. 2:74 to 1:5, as previously demonstrated in the literature for true hard spheres. A number of simple analytic expressions that represent quite well the percolation characteristics of the WCA system are proposed.

Robust and accurate method for free-energy calculation of charged molecular systems.

A new approach is presented to eliminate the problem of creation and/or annihilation of atoms in free-energy calculations of charged molecular systems. The method employs a damping potential in the Ewald summation scheme, which is an exact solution of the electrostatics for three-dimensional periodic systems. The proposed method enables entire molecules to be mutated from a noninteracting (ideal) state in an efficient and robust way, thus providing a means by which accurate absolute free energies of structurally complex molecules can be determined. This methodology will enable chemical and phase equilibria to be determined for large molecular species with significant charge distributions, e.g., biomolecules and drugs.