Melting line of the Lennard-Jones system, infinite size, and full potential

Literature estimates of the melting curve of the Lennard-Jones system vary by as much as 10%. The origin of such discrepancies remains unclear. We present precise values for the Lennard-Jones melting temperature, and we examine possible sources of systematic errors in the prediction of melting points, including finite-size and interaction-cutoff effects. A hypothetical thermodynamic integration path is used to find the relative free energies of the solid and liquid phases, for various system sizes, at constant cutoff radius. The solid-liquid relative free energy and melting temperature scale linearly as the inverse of the number of particles, and it is shown that finite-size effects can account for deviations in the melting temperature (from the infinite-size limit) of up to 5%. An extended-ensemble density-of-states method is used to determine free energy changes for each phase as a continuous function of the cutoff radius. The resulting melting temperature predictions exhibit an oscillatory behavior as the cutoff radius is increased. Deviations in the melting temperature (from the full potential limit) arising from a finite cutoff radius are shown to be of comparable magnitude as those resulting from finite-size effects. This method is used to identify melting temperatures at five different pressures, for the infinite-size and full potential Lennard-Jones system. We use our simulation results as references to connect the Lennard-Jones solid equation of state of van der Hoef with the Lennard-Jones fluid equation of state of Johnson. Once the references are applied the two equations of state are used to identify a melting curve. An empirical equation that fits this melting curve is provided. We also report a reduced triple point temperature T(tr)=0.694.     

Relationships between self-diffusivity, packing fraction, and excess entropy in simple bulk and confined fluids

Static measures such as density and entropy, which are intimately connected to structure, have featured prominently in modern thinking about the dynamics of the liquid state. Here, we explore the connections between self-diffusivity, density, and excess entropy for two of the most widely used model "simple" liquids, the equilibrium Lennard-Jones and square-well fluids, in both bulk and confined environments. We find that the self-diffusivity data of the Lennard-Jones fluid can be approximately collapsed onto a single curve (i) versus effective packing fraction and (ii) in appropriately reduced form versus excess entropy, as suggested by two well-known scaling laws. Similar data collapse does not occur for the square-well fluid, a fact that can be understood on the basis of the nontrivial effects that temperature has on its static structure. Nonetheless, we show that the implications of confinement for the self-diffusivity of both of these model fluids, over a broad range of equilibrium conditions, can be predicted on the basis of knowledge of the bulk fluid behavior and either the effective packing fraction or the excess entropy of the confined fluid. Excess entropy is perhaps the most preferable route due to its superior predictive ability and because it is a standard, unambiguous thermodynamic quantity that can be readily predicted via classical density functional theories of inhomogeneous fluids.     

Energy conservation in molecular dynamics simulations of classical systems

Classical Newtonian dynamics is analytic and the energy of an isolated system is conserved. The energy of such a system, obtained by the discrete "Verlet" algorithm commonly used in molecular dynamics simulations, fluctuates but is conserved in the mean. This is explained by the existence of a "shadow Hamiltonian" H [S. Toxvaerd, Phys. Rev. E 50, 2271 (1994)], i.e., a Hamiltonian close to the original H with the property that the discrete positions of the Verlet algorithm for H lie on the analytic trajectories of H. The shadow Hamiltonian can be obtained from H by an asymptotic expansion in the time step length. Here we use the first non-trivial term in this expansion to obtain an improved estimate of the discrete values of the energy. The investigation is performed for a representative system with Lennard-Jones pair interactions. The simulations show that inclusion of this term reduces the standard deviation of the energy fluctuations by a factor of 100 for typical values of the time step length. Simulations further show that the energy is conserved for at least one hundred million time steps provided the potential and its first four derivatives are continuous at the cutoff. Finally, we show analytically as well as numerically that energy conservation is not sensitive to round-off errors.     

Canonical sampling through velocity rescaling

The authors present a new molecular dynamics algorithm for sampling the canonical distribution. In this approach the velocities of all the particles are rescaled by a properly chosen random factor. The algorithm is formally justified and it is shown that, in spite of its stochastic nature, a quantity can still be defined that remains constant during the evolution. In numerical applications this quantity can be used to measure the accuracy of the sampling. The authors illustrate the properties of this new method on Lennard-Jones and TIP4P water models in the solid and liquid phases. Its performance is excellent and largely independent of the thermostat parameter also with regard to the dynamic properties.     

Long-range contribution to inverse-distance power and Morse energy and pressure terms of a molecular liquid,tested on benzene

Analytical expressions were derived to estimate the contribution to the intermolecular energy and pressure from pairs of rigid molecules farther than a cutoff distance; inverse power terms (like in the Lennard-Jones potential) and exponential terms (like in the Morse potential) were considered. The Lennard-Jones case was tested on liquid benzene at room temperature in the Gibbs ensemble using the cavity-biased technique. The density produced by simulations using the derived cutoff correction was found to be independent of the cutoff used.     

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.     

Direct calculation of solid-liquid coexistence points of a binary mixture by thermodynamic integration

We present a new thermodynamic integration method that directly connects the liquid and the solid phases of a binary mixture by a reversible path. The states along the path are simulated in the isothermal-isobaric semigrand canonical ensemble, in which temperature, pressure, the total number of particles, and the fugacity fractions of the components are held fixed. The thermodynamic integration yields the chemical-potential difference between the two phases for one of the components and this information is then used to locate the solid-liquid coexistence points. The melting temperatures predicted by our method agree well with those predicted by the Gibbs-Duhem integration for a truncated and shifted Lennard-Jones system with a cutoff radius of 2.5sigma.     

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.     

Beyond the Lennard-Jones model: a simple and accurate potential function probed by high resolution scattering data useful for molecular dynamics simulations

Scattering data, measured for rare gas-rare gas systems under high angular and energy resolution conditions, have been used to probe the reliability of a recently proposed interaction potential function, which involves only one additional parameter with respect to the venerable Lennard-Jones (LJ) model and is hence called Improved Lennard-Jones (ILJ). The ILJ potential eliminates most of the inadequacies at short- and long-range of the LJ model. Further reliability tests have been performed by comparing calculated vibrational spacings with experimental values and calculated interaction energies at short-range with those obtained from the inversion of gaseous transport properties. The analysis, extended also to systems involving ions, suggests that the ILJ potential model can be used to estimate the behavior of unknown systems and can help to assess the different role of the leading interaction components. Moreover, due to its simple formulation, the physically reliable ILJ model appears to be particularly useful for molecular dynamics simulations of both neutral and ionic systems.     

An improved algorithm to locate critical points in a 3D scalar field as implemented in the program MORPHY

A new algorithm for location of the critical points in general scalar fields is described. The new method has been developed as part of an on-going process to exploit the topologic analysis of general 3D scalar fields. Part of this process involves the use of topologic information to seed the critical point search algorithm. The continuing move away from topologic studies of just the electron density requires more general algorithms and the ability to easily "plug in" new functions, for example, the Laplacian of the electron density ( triangle down (2)rho), the Electron Localisation Function (ELF), the Localised Orbital Locator (LOL), the Lennard-Jones function (LJF), as well as any new functions that may be proposed in the future. Another important aspect of the current algorithm is the retention of all possible intermediate information, for example, the paths describing the connectivity of critical points, as well as an ability to restart searches, something that becomes increasingly important when analysing larger systems. This new algorithm represents a core part of a new local version of the MORPHY code. We distinguish nine universal types of gradient paths.     

Toward a robust and general molecular simulation method for computing solid-liquid coexistence

A rigorous and generally applicable method for computing solid-liquid coexistence is presented. The method overcomes some of the technical difficulties associated with other solid-liquid simulation procedures and can be implemented within either a molecular dynamics or Monte Carlo framework. The method consists of three steps: First, relative Gibbs free energy curves are created for the solid and liquid phases using histogram reweighting. Next, the free energy difference between the solid and liquid phases is evaluated at a single state point by integrating along a pseudosupercritical transformation path that connects the two phases. Using this result, the solid and liquid free energy curves are referenced to a common point, allowing a single coexistence point to be determined. Finally, Gibbs-Duhem integration is used to determine the full coexistence curve. To evaluate its utility, this method is applied to the Lennard-Jones and NaCl systems. Results for solid-liquid coexistence agree with previous calculations for these systems. In addition, it is shown that the NaCl model does not correctly describe solid-liquid coexistence at high pressures. An analysis of the accuracy of the method indicates that the results are most sensitive to the transformation free energy calculation.     

Global minimum structure searches via particle swarm optimization

Novel implementation of the evolutionary approach known as particle swarm optimization (PSO) capable of finding the global minimum of the potential energy surface of atomic assemblies is reported. This is the first time the PSO technique has been used to perform global optimization of minimum structure search for chemical systems. Significant improvements have been introduced to the original PSO algorithm to increase its efficiency and reliability and adapt it to chemical systems. The developed software has successfully found the lowest-energy structures of the LJ(26) Lennard-Jones cluster, anionic silicon hydride Si(2)H(5) (-), and triply hydrated hydroxide ion OH(-) (H(2)O)(3). It requires relatively small population sizes and demonstrates fast convergence. Efficiency of PSO has been compared with simulated annealing, and the gradient embedded genetic algorithm.     

Application of Ewald summations to long-range dispersion forces

Results illustrating the effects of using explicit summation terms for the r(-6) dispersion term on the interfacial properties of a Lennard-Jones fluid and SPC/E water are presented. For the Lennard-Jones fluid, we find that the use of long-range summations, even with a short "crossover radius," yields results that are consistent with simulations using large cutoff radii. Simulations of SPC/E water demonstrate that the long-range dispersion forces are of secondary importance to the Coulombic forces. In both cases, we find that the ratio of the box size L( parallel) to the crossover radius r(c) (k) plays an important role in determining the magnitude of the long-range dispersion correction, although its effect is secondary when Coulombic interactions are also present.     

Communication: Shifted forces in molecular dynamics

Simulations involving the Lennard-Jones potential usually employ a cutoff at r = 2.5σ. This communication investigates the possibility of reducing the cutoff. Two different cutoff implementations are compared, the standard shifted potential cutoff and the less commonly used shifted forces cutoff. The first has correct forces below the cutoff, whereas the shifted forces cutoff modifies Newton's equations at all distances. The latter is nevertheless superior; we find that for most purposes realistic simulations may be obtained using a shifted forces cutoff at r = 1.5σ, even though the pair force is here 30 times larger than at r = 2.5σ.     

Is there a relation between excess volume and miscibility in binary liquid mixtures?

Molecular dynamics computer simulations of various symmetrical Lennard-Jones (LJ) models are used to elucidate how the excess volume in dense binary liquids is related to the microscopic interactions between the particles. Both fully miscible systems and systems with a liquid-liquid phase separation are considered by varying systematically the parameters of the LJ potentials. The phase diagrams with the critical points of the demixing systems are determined by means of Monte Carlo simulations in the semigrandcanonical ensemble. The different LJ models are investigated by computing Bhatia-Thornton structure factors, enthalpy of mixing, and excess volume. For the demixing systems, the LJ models show a positive enthalpy of mixing while it is negative for the systems without miscibility gap. In contrast to that, the excess volume can be negative and positive for both demixing and fully miscible systems. This behavior is explained in terms of the interplay between the repulsive and attractive terms in the LJ potential. Whereas repulsions dominate the packing of particles as reflected by the number-density structure factor, the chemical ordering and thus the concentration structure factor are strongly affected by attractive interactions, leading to the "anomalies" of the excess volume.     

Stochastic potential switching algorithm for Monte Carlo simulations of complex systems

This paper describes a new Monte Carlo method based on a novel stochastic potential switching algorithm. This algorithm enables the equilibrium properties of a system with potential V to be computed using a Monte Carlo simulation for a system with a possibly less complex stochastically altered potential V. By proper choices of the stochastic switching and transition probabilities, it is shown that detailed balance can be strictly maintained with respect to the original potential V. The validity of the method is illustrated with a simple one-dimensional example. The method is then generalized to multidimensional systems with any additive potential, providing a framework for the design of more efficient algorithms to simulate complex systems. A near-critical Lennard-Jones fluid with more than 20,000 particles is used to illustrate the method. The new algorithm produced a much smaller dynamic scaling exponent compared to the Metropolis method and improved sampling efficiency by over an order of magnitude.     

Direct calculation of solid-vapor coexistence points by thermodynamic integration: application to single component and binary systems

We present a new thermodynamic integration method that directly connects the vapor and solid phases by a reversible path. The thermodynamic integration in the isothermal-isobaric ensemble yields the Gibbs free energy difference between the two phases, from which the sublimation temperature can be easily calculated. The method extends to the binary mixture without any modification to the integration path simply by employing the isothermal-isobaric semigrand ensemble. The thermodynamic integration, in this case, yields the chemical potential difference between the solid and vapor phases for one of the components, from which the binary sublimation temperature can be calculated. The coexistence temperatures predicted by our method agree well with those in the literature for single component and binary Lennard-Jones systems.     

Communication: phase transitions, criticality, and three-phase coexistence in constrained cell models

In simulation studies of fluid-solid transitions, the solid phase is usually modeled as a constrained system in which each particle is confined to move in a single Wigner-Seitz cell. The constrained cell model has been used in the determination of fluid-solid coexistence via thermodynamic integration and other techniques. In the present work, the phase diagram of such a constrained system of Lennard-Jones particles is determined from constant-pressure simulations. The pressure-density isotherms exhibit inflection points which are interpreted as the mechanical stability limit of the solid phase. The phase diagram of the constrained system contains a critical and a triple point. The temperature and pressure at the critical and the triple point are both higher than those of the unconstrained system due to the reduction in the entropy caused by the single occupancy constraint.