NVU dynamics. I. Geodesic motion on the constant-potential-energy hypersurface

An algorithm is derived for computer simulation of geodesics on the constant-potential-energy hypersurface of a system of N classical particles. First, a basic time-reversible geodesic algorithm is derived by discretizing the geodesic stationarity condition and implementing the constant-potential-energy constraint via standard Lagrangian multipliers. The basic NVU algorithm is tested by single-precision computer simulations of the Lennard-Jones liquid. Excellent numerical stability is obtained if the force cutoff is smoothed and the two initial configurations have identical potential energy within machine precision. Nevertheless, just as for NVE algorithms, stabilizers are needed for very long runs in order to compensate for the accumulation of numerical errors that eventually lead to "entropic drift" of the potential energy towards higher values. A modification of the basic NVU algorithm is introduced that ensures potential-energy and step-length conservation; center-of-mass drift is also eliminated. Analytical arguments confirmed by simulations demonstrate that the modified NVU algorithm is absolutely stable. Finally, we present simulations showing that the NVU algorithm and the standard leap-frog NVE algorithm have identical radial distribution functions for the Lennard-Jones liquid.     

Dynamical properties of the soft sticky dipole-quadrupole-octupole water model: a molecular dynamics study

The dynamical properties of the soft sticky dipole-quadrupole-octupole (SSDQO) water model using SPC/E moments are calculated utilizing molecular dynamics simulations. This new potential for liquid water describes the water-water interactions by a Lennard-Jones term and a sticky potential, which is an approximate moment expansion with point dipole, quadrupole, and octupole moments, and reproduces radial distribution functions of pure liquid water using the moments of SPC/E [Ichiye and Tan, J. Chem. Phys. 124, 134504 (2006)]. The forces and torques of SSDQO water for the dipole-quadrupole, quadrupole-quadrupole, and dipole-octupole interactions are derived here. The simulations are carried out at 298 K in the microcanonical ensemble employing the Ewald method for the long-range dipole-dipole interactions. Here, various dynamical properties associated with translational and rotational motions of SSDQO water using the moments of SPC/E (SSDQO:SPC/E) water are compared with the results from SPC/E and also experiment. The self-diffusion coefficient of SSDQO:SPC/E water is found to be in excellent agreement with both SPC/E and experiment whereas the single particle orientational relaxation time for dipole vector is better than SPC/E water but it is somewhat smaller than experiment. The dielectric constant of SSDQO:SPC/E is essentially identical to SPC/E, and both are slightly lower than experiment. Also, molecular dynamics simulations of the SSDQO water model are found to be about twice as fast as three-site models such as SPC/E.     

On combining molecular dynamics and stochastic dynamics simulations to compute reaction rates in liquids

An approach that combines molecular dynamics and stochastic dynamics calculations for obtaining reaction rates in liquids is investigated by studying the cis-->trans isomerization of HONO in liquid krypton. The isomerization rates are computed for several liquid densities by employing full-dimensional molecular-dynamics simulations. The rates are also computed by employing the stochastic dynamics method for a wide range of collision frequencies. Comparisons of the two sets of the computed rates show that for a wide range of liquid densities there is a simple linear relation between the liquid density rho and the collision frequency alpha, that is, alpha=crho. This suggests that once the constant c is determined from a molecular-dynamics calculation at a single density, the reaction rates can be obtained from stochastic dynamics calculations for the entire range of liquid densities where alpha=crho holds. The applicability of the combined molecular dynamics and stochastic dynamics approach provides a practical means for obtaining rate constants at considerable savings of computer time compared to that required by using full-dimensional molecular-dynamics simulations alone.     

Molecular dynamics study of nanoparticle stability at liquid interfaces: effect of nanoparticle-solvent interaction and capillary waves

While the interaction of colloidal particles (sizes in excess of 100 nm) with liquid interfaces may be understood in terms of continuum models, which are grounded in macroscopic properties such as surface and line tensions, the behaviour of nanoparticles at liquid interfaces may be more complex. Recent simulations [D. L. Cheung and S. A. F. Bon, Phys. Rev. Lett. 102, 066103 (2009)] of nanoparticles at an idealised liquid-liquid interface showed that the nanoparticle-interface interaction range was larger than expected due, in part, to the action of thermal capillary waves. In this paper, molecular dynamics simulations of a Lennard-Jones nanoparticle in a binary Lennard-Jones mixture are used to confirm that these previous results hold for more realistic models. Furthermore by including attractive interactions between the nanoparticle and the solvent, it is found that the detachment energy decreases as the nanoparticle-solvent attraction increases. Comparison between the simulation results and recent theoretical predictions [H. Lehle and M. Oettel, J. Phys. Condens. Matter 20, 404224 (2008)] shows that for small particles the incorporation of capillary waves into the predicted effective nanoparticle-interface interaction improves agreement between simulation and theory.     

Cavitation and crystallization in a metastable Lennard-Jones liquid at negative pressures and low temperatures

Molecular dynamics simulations have been used to investigate the kinetics of spontaneous cavitation and crystallization in a Lennard-Jones liquid at negative pressures in the temperature range where these processes compete with each other. The nucleation rate has been calculated in NVE and NpT ensembles by the method of mean lifetime and the transition interface sampling method with parallel path swapping. The data obtained have been used to determine in the framework of classical nucleation theory the value of the ratio of the solid-liquid and the liquid-void interfacial free energy for critical crystals and cavities and the values of their volumes at points where the cavitation rate of the liquid is equal to the rate of its crystallization.     

Role of the first coordination shell in determining the equilibrium structure and dynamics of simple liquids

The traditional view that the physical properties of a simple liquid are determined primarily by its repulsive forces was recently challenged by Berthier and Tarjus, who showed that in some cases ignoring the attractions leads to large errors in the dynamics [L. Berthier and G. Tarjus, Phys. Rev. Lett. 103, 170601 (2009); J. Chem. Phys. 134, 214503 (2011)]. We present simulations of the standard Lennard-Jones liquid at several condensed-fluid state points, including a fairly low density state and a very high density state, as well as simulations of the Kob-Andersen binary Lennard-Jones mixture. By varying the range of the forces via a shifted-forces cutoff, results for the thermodynamics, dynamics, and structure show that the determining factor for getting the correct statics and dynamics is not whether or not the attractive forces per se are included in the simulations. What matters is whether or not interactions are included from all particles within the first coordination shell - the attractive forces can thus be ignored, but only at extremely high densities. The recognition of the importance of a local shell in condensed fluids goes back to van der Waals; our results confirm this idea and thereby the basic picture of the old hole and cell theories for simple condensed fluids.     

Tests of an approximate scaling principle for dynamics of classical fluids

We used molecular dynamics computer simulations to test an approximate scaling principle that conjectures that two equilibrium atomic liquids have very similar dynamical properties if they have the same density and similar static pair correlation functions when the length scales of the two liquids are adjusted appropriately, even if they have different interatomic potentials and different temperatures. The simulations were performed on two types of model atomic liquids at various temperatures at the same density. In the first type, the interatomic potential is the Lennard-Jones potential (LJ). In the second type, the interatomic potential is the repulsive part of the Lennard-Jones potential (RLJ). We identified pairs of systems that have very similar pair correlation functions despite the fact that they had different potentials. Each pair consisted of an LJ liquid at a specific temperature and a corresponding RLJ liquid at a lower temperature. We compared various time correlation functions and transport coefficients of the two systems in each pair. Many dynamical properties are very similar in each pair, in accordance with the approximate scaling principle, whereas others are significantly different. The results indicate that certain dynamical properties are very insensitive to large changes in the interatomic potential that leave the pair correlation function largely unchanged, whereas other dynamical properties are much more sensitive to such changes in the potential. The transport coefficients for diffusion and viscosity are among the dynamical properties that are insensitive to such changes in the potential, and this may be part of the reason transport properties of many fluids have been calculated or rationalized in terms of a simple hard sphere model of liquids.     

Discontinuous molecular dynamics for rigid bodies: applications

Event-driven molecular dynamics simulations are carried out on two rigid-body systems which differ in the symmetry of their molecular mass distributions. First, simulations of methane in which the molecules interact via discontinuous potentials are compared with simulations in which the molecules interact through standard continuous Lennard-Jones potentials. It is shown that under similar conditions of temperature and pressure, the rigid discontinuous molecular dynamics method reproduces the essential dynamical and structural features found in continuous-potential simulations at both gas and liquid densities. Moreover, the discontinuous molecular dynamics approach is demonstrated to be between 3 and 100 times more efficient than the standard molecular dynamics method depending on the specific conditions of the simulation. The rigid discontinuous molecular dynamics method is also applied to a discontinuous-potential model of a liquid composed of rigid benzene molecules, and equilibrium and dynamical properties are shown to be in qualitative agreement with more detailed continuous-potential models of benzene. The few qualitative differences in the angular dynamics of the two models are related to the relatively crude treatment of variations in the discontinuous repulsive interactions as one benzene molecule rotates by another.     

Structure and diffusion of nanoparticle monolayers floating at liquid/vapor interfaces: a molecular dynamics study

Large-scale molecular dynamics simulations are used to simulate a layer of nanoparticles floating on the surface of a liquid. Both a low viscosity liquid, represented by Lennard-Jones monomers, and a high viscosity liquid, represented by linear homopolymers, are studied. The organization and diffusion of the nanoparticles are analyzed as the nanoparticle density and the contact angle between the nanoparticles and liquid are varied. When the interaction between the nanoparticles and liquid is reduced the contact angle increases and the nanoparticles ride higher on the liquid surface, which enables them to diffuse faster. In this case the short-range order is also reduced as seen in the pair correlation function. For the polymeric liquids, the out-of-layer fluctuation is suppressed and the short-range order is slightly enhanced. However, the diffusion becomes much slower and the mean square displacement even shows sub-linear time dependence at large times. The relation between diffusion coefficient and viscosity is found to deviate from that in bulk diffusion. Results are compared to simulations of the identical nanoparticles in 2-dimensions.     

Molecular dynamics simulation of liquid N2O4<-->2NO2 by orientation-sensitive pairwise potential. III. Reaction dynamics

The dissociation and association dynamics of N2O4 [see text] 2NO2 in liquid state are studied by classical molecular dynamics simulations of reactive liquid NO2. An OSPP+LJ potential between NO2 molecules, which is a sum of an orientation-sensitive pairwise potential (OSPP) between N-N atoms proposed in Paper I [J. Chem. Phys. 115, 10852 (2001)] and Lennard-Jones potentials between N-O and O-O atoms, has been used in the simulation. The reaction dynamics is studied as a function of well depth De and anisotropy factors of the OSPP potential: Atheta (0< or =Atheta< or =1) for the rocking angle and Atau (0< or =Atau< or =0.5) for the torsional angle of relative NO2-NO2 orientation. The lifetime tauD of initially prepared NO2 dimers is found to increase as De increases, Atheta increases, and Atau decreases. Dissociation and association dynamics are studied in detail around the extreme limit of pure NO2-dimer liquid: De=0.12 x 10(-18) J, Atheta=0.5, and Atau=0.1, which has been found to reproduce both the observed liquid phase equilibrium properties and Raman band shapes of the dissociation mode very well. The dissociation dynamics from microscopic reaction trajectories is compared with the potential of the mean force (PMF) as a function of the N-N distance R. The PMF of reactive liquid NO2 shows a transition state barrier at R=2.3-2.5 A, and NO2-trimer structure is found to be formed at the barrier. Two types of dissociation of the NO2 dimer-the dissociation by collisional activation of the reactive mode to cross the dissociation limit and the NO2-mediated dissociation via bond transfer-are studied. The latter needs less free energy and is found to be much more probable. The dissociation trajectories and PMF in reactive liquid NO2 are compared with those of a reactive NO2 pair in inert solvent N2O4.     

Molecular dynamics simulation of potassium along the liquid-vapor coexistence curve

The applicability of pair potential functions to liquid alkali metals is questionable. On the one hand, some recent reports in the literature suggest the validity of two-parameter pair-wise additive Lennard-Jones (LJ) potentials for liquid alkali metals. On the other hand, there are some reports suggesting the inaccuracy of pair potential functions for liquid metals. In this work, we have performed extensive molecular dynamics simulations of vapor-liquid phase equilibria in potassium to check the validity of the proposed LJ potentials and to improve their accuracy by changing the LJ exponents and taking into account the temperaturedependencies of the potential parameters. We have calculated the orthobaric liquid and vapor densities of potassium using LJ (12–6), LJ (8.5–4) and LJ (5–4), effective pair potential energy functions. The results show that using an LJ (8.5–4) potential energy function with temperature-independent parameters, ε and σ, is inadequate to account for the vapor-liquid coexistence properties of potassium. Taking into account the temperature-dependencies of the LJ parameters, ε(T) and σ(T), we obtained the densities of coexisting liquid and vapor potassium in a much better agreement with experimental data. Changing the magnitude of repulsive and attractive contributions to the potential energy function shows that a two-parameter LJ (5–4) potential can well reproduce the densities of liquid and vapor potassium. The results show that LJ (5–4) potential with temperature-dependent parameters produces the densities of liquid and vapor potassium more accurately, compared to the results obtained using LJ (12–6) and LJ (8.5–4) potential energy functions.     

Molecular dynamics simulation of the energetic room-temperature ionic liquid, 1-hydroxyethyl-4-amino-1,2,4-triazolium nitrate (HEATN)

Molecular dynamics (MD) simulations have been performed to investigate the structure and dynamics of an energetic ionic liquid, 1-hydroxyethyl-4-amino-1,2,4-triazolium nitrate (HEATN). The generalized amber force field (GAFF) was used, and an electronically polarizable model was further developed in the spirit of our previous work (Yan, T.; Burnham, C. J.; Del Popolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 11877). In the process of simulated annealing from a liquid state at 475 K down to a glassy state at 175 K, the MD simulations identify a glass-transition temperature region at around 250-275 K, in agreement with experiment. The self-intermediate scattering functions show vanishing boson peaks in the supercooled region, indicating that HEATN may be a fragile glass former. The coupling/decoupling of translational and reorientational ion motion is also discussed, and various other physical properties of the liquid state are intensively studied at 400 K. A complex hydrogen bond network was revealed with the calculation of partial radial distribution functions. When compared to the similarly sized 1-ethyl-4-methyl-1,4-imidazolium nitrate ionic liquid, EMIM+/NO3-, a hydrogen bond network directly resulting in the poorer packing efficiency of ions is observed, which is responsible for the lower melting/glass-transition point. The structural properties of the liquid/vacuum interface shows that there is vanishing layering at the interface, in accordance with the poor ion packing. The effects of electronic polarization on the self-diffusion, viscosity, and surface tension of HEATN are found to be significant, in agreement with an earlier study on EMIM+/NO3- (Yan, T.; Burnham, C. J.; Del Popolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 11877).     

Pressure-energy correlations in liquids. V. Isomorphs in generalized Lennard-Jones systems

This series of papers is devoted to identifying and explaining the properties of strongly correlating liquids, i.e., liquids with more than 90% correlation between their virial W and potential energy U fluctuations in the NVT ensemble. Paper IV [N. Gnan et al., J. Chem. Phys. 131, 234504 (2009)] showed that strongly correlating liquids have "isomorphs," which are curves in the phase diagram along which structure, dynamics, and some thermodynamic properties are invariant in reduced units. In the present paper, using the fact that reduced-unit radial distribution functions are isomorph invariant, we derive an expression for the shapes of isomorphs in the WU phase diagram of generalized Lennard-Jones systems of one or more types of particles. The isomorph shape depends only on the Lennard-Jones exponents; thus all isomorphs of standard Lennard-Jones systems (with exponents 12 and 6) can be scaled onto a single curve. Two applications are given. One tests the prediction that the solid-liquid coexistence curve follows an isomorph by comparing to recent simulations by Ahmed and Sadus [J. Chem. Phys. 131, 174504 (2009)]. Excellent agreement is found on the liquid side of the coexistence curve, whereas the agreement is less convincing on the solid side. A second application is the derivation of an approximate equation of state for generalized Lennard-Jones systems by combining the isomorph theory with the Rosenfeld-Tarazona expression for the temperature dependence of the potential energy on isochores. It is shown that the new equation of state agrees well with simulations.     

A box-counting-based algorithm for computing Shannon entropy in molecular dynamics simulations

A box-counting-based algorithm (SEBC) has been developed for the numerical computation of the Shannon entropy from samples of continuous functions. Its performance was tested by applying it to several samples of known continuous distribution functions. The results obtained with SEBC reproduced those obtained by analytical or numerical integration. SEBC was also employed for computing the Shannon entropies of the steric energy, Sh(E(S)), of several amino acids from their in vacuo NVE molecular dynamics simulations using the AMBER-4 force field. The results obtained correlate linearly with the experimental standard thermodynamic entropies of these compounds. This work points to the possibility of introducing straightforward and reliable calculations of thermodynamic entropies from empirical linear relationships with Sh(E(S)) obtained from MD simulations.     

A molecular dynamics study of nitric oxide in water: diffusion and structure

We present molecular dynamics simulations of the diffusion coefficients and structure of water-nitric oxide mixtures at ambient (298 K) and in vivo (310 K) conditions. A two-site rigid-body molecular model with partial charges and a Lennard-Jones potential on both sites is proposed for nitric oxide and used in conjunction with the extended simple point-charge model for liquid water in our simulations. The diffusion coefficients obtained from the simulations are in good agreement with experimental data. The results from intermolecular partial pair functions show that under these thermodynamic conditions, the existence of nitric oxide in liquid water has little impact on the structure of water and the tendency to form H bonds between water molecules. We also find that it is unlikely that H bonds form between the hydrogen atoms in water and either the nitrogen or the oxygen atom on the nitric oxide at the temperatures and densities examined in this study. This study suggests that in low concentrations nitric oxide molecules exist as free molecules in liquid water rather than forming complexes with water molecules.     

A molecular dynamics study of water nucleation using the TIP4P/2005 model

Extensive molecular dynamics simulations were conducted using the TIP4P/2005 water model of Abascal and Vega [J. Chem. Phys. 123, 234505 (2005)] to investigate its condensation from supersaturated vapor to liquid at 330 K. The mean first passage time method [J. Wedekind, R. Strey, and D. Reguera, J. Chem. Phys. 126, 134103 (2007); L. S. Bartell and D. T. Wu, 125, 194503 (2006)] was used to analyze the influence of finite size effects, thermostats, and charged species on the nucleation dynamics. We find that the Nose?-Hoover thermostat and the one proposed by Bussi et al. [J. Chem. Phys. 126, 014101 (2007)] give essentially the same averages. We identify the maximum thermostat coupling time to guarantee proper thermostating for these simulations. The presence of charged species has a dramatic impact on the dynamics, inducing a marked change towards a pure growth regime, which highlights the importance of ions in the formation of liquid droplets in the atmosphere. It was found a small but noticeable sign preference at intermediate cluster sizes (between 5 and 30 water molecules) corresponding mostly to the formation of the second solvation shell around the ion. The TIP4P/2005 water model predicts that anions induce faster formation of water clusters than cations of the same magnitude of charge.     

Simulating critical dynamics in liquid mixtures: short-range and long-range contributions

Recently, Das et al. [J. Chem. Phys. 125, 024506 (2006)] established that computer simulations of critical dynamics in a binary Lennard-Jones mixture are consistent with the predicted Stokes-Einstein behavior of the asymptotic decay rate of the order-parameter fluctuations near criticality. Here, we show that the noncritical or "background" contributions to the computed diffusion coefficient are also in agreement with both theory and experiment, thus further validating the feasibility of molecular dynamics simulations for studying dynamic critical behavior.     

Pair correlation functions in mixtures of Lennard-Jones particles

The pair correlation functions for a mixture of two Lennard-Jones particles were computed by both the Percus-Yevick equations and by molecular dynamics. The changes in the pair correlation function resulting from changes in the composition of the mixtures are quite unexpected. Essentially, identical changes are obtained from the Percus-Yevick equations and from molecular dynamics simulations. The molecular reason for this unexpected behavior is discussed.     

The dynamics of ions in dimethylformamide

Time autocorrelation functions were calculated from the results of molecular dynamics simulations of ions in liquid N,N-dimethylformamide at 233, 298, and 318 K. The mobilities and diffusion coefficients of ions determined by molecular dynamics simulations and from conductometric data are compared. The influence of the temperature and size of ions on the mobility of solvent molecules in the first solvation shells of ions is analyzed.     

Spontaneous and induced dynamic fluctuations in glass formers. I. General results and dependence on ensemble and dynamics

We study theoretically and numerically a family of multipoint dynamic susceptibilities that quantify the strength and characteristic length scales of dynamic heterogeneities in glass-forming materials. We use general theoretical arguments (fluctuation-dissipation relations and symmetries of relevant dynamical field theories) to relate the sensitivity of averaged two-time correlators to temperature and density to spontaneous fluctuations of the local dynamics. Our theoretical results are then compared to molecular dynamics simulations of the Newtonian, Brownian, and Monte Carlo dynamics of two representative glass-forming liquids, a fragile binary Lennard-Jones mixture, and a model for the strong glass-former silica. We justify in detail the claim made by Berthier et al. [Science 310, 1797 (2005)] that the temperature dependence of correlation functions allows one to extract useful information on dynamic length scales in glassy systems. We also discuss some subtle issues associated with the choice of microscopic dynamics and of statistical ensemble through conserved quantities, which are found to play an important role in determining dynamic correlations.