A classical polarizable model for simulations of water and ice

We develop a classical rigid polarizable model of water for molecular simulations of water and ice. The model uses the Rowlinson five-site geometry: oxygen bearing the Lennard-Jones interaction and linearly polarizable point dipole, two positively charged hydrogens, and two massless negative charges placed symmetrically off oxygen so that the experimental dipole moment is reproduced. The target properties are the densities of water and ice Ih, diffusivity, enthalpies of fusion and vaporization, and the ice Ih melting point. The surface tension at lower temperatures is by 7% underestimated whereas the dielectric constant by 6% overestimated. Diffusivity and viscosity worsen at higher temperatures, although the Stokes radius is overestimated only by 2-7%. The ice Ih melting temperature is 260 K and the temperature of maximum density is 269 K. Rescaling the charges by a factor of 1.01 and Lennard-Jones energy by 1.0201 improves the melting point and energy-related quantities but shifts the agreement of kinetic properties to higher temperatures. For the model we propose abbreviation POL4D.     

Molecular dynamics simulations of vapor/liquid coexistence using the nonpolarizable water models

The surface tension, vapor-liquid equilibrium densities, and equilibrium pressure for common water models were calculated using molecular dynamics simulations over temperatures ranging from the melting to the critical points. The TIP4P/2005 and TIP4P-i models produced better values for the surface tension than the other water models. We also examined the correlation of the data to scaling temperatures based on the critical and melting temperatures. The reduced temperature (T/T(c)) gives consistent equilibrium densities and pressure, and the shifted temperature T + (T(c, exp) - T(c, sim)) gives consistent surface tension among all models considered in this study. The modified fixed charge model which has the same Lennard-Jones parameters as the TIP4P-FQ model but uses an adjustable molecular dipole moment is also simulated to find the differences in the vapor-liquid coexistence properties between fixed and fluctuating charge models. The TIP4P-FQ model (2.72 Debye) gives the best estimate of the experimental surface tension. The equilibrium vapor density and pressure are unaffected by changes in the dipole moment as well as the surface tension and liquid density.     

Thermal fluctuations of clusters with the long-range interaction

Analysis of surface fluctuation spectra is performed for a large cluster of particles interacting via a sum of the short-range Lennard-Jones potential and long-range ±1/r potential, where the positive sign corresponds to the gravity, and negative corresponds to the electrostatic interaction. The spectral amplitudes of thermally driven capillary modes in a self-consistent field induced by cluster particles including the modes with no axial symmetry are derived in the approximation of small amplitudes. It is demonstrated that within used approximation, the surface tension is independent of the field strength. The low wave vector amplitudes are damped by attracting field that compresses the cluster and magnified by repulsing field leading to cluster fission. The fission threshold is found to be different from that found by Bohr and Wheeler and Frenkel due to the replacement of the ordinary surface tension by the bare one. Molecular dynamics study of a cluster with the long-range interaction in the vapor environment is performed using a novel integrator for a multiscale system. Simulation scheme implies rotation of the long-range components of forces acting on cluster particles thus vanishing an artificial torque. Simulation results justify theoretical conclusion of modes damping and independence of the surface tension of the field strength. Fission threshold evaluated from simulation data is in a good agreement with theory.     

Surface tension of water droplets upon homogeneous droplet nucleation in water vapor

A method has been proposed for determining interfacial free energy from the data of molecular dynamics simulation. The method is based on the thermodynamic integration procedure and is distinguished by applicability to both planar interfaces and those characterized by a high curvature. The workability of the method has been demonstrated by the example of determining the surface tension for critical nuclei of water droplets upon condensation of water vapor. The calculation has been performed at temperatures of 273–373 K and a pressure of 1 atm, thus making it possible to determine the temperature dependence of the surface tension for water droplets and compare the results obtained with experimental data and the simulation results for a “planar” vapor–liquid interface.     

On the long-range corrections to computer simulation results for the Lennard-Jones vapor-liquid interface

The long-range corrections (LRCs) to the configurational energy have been taken into consideration in the Monte Carlo simulation of the vapor-liquid interface for a pure Lennard-Jones (LJ) fluid. The simulated bulk densities agree satisfactorily with those obtained from the Gibbs ensemble method, and the simulated surface tension values agree reasonably well with those reported in the literature for a larger number of molecules and a larger cut-off distance. To compare the influence of the potential forms on the simulation results, a truncated LJ potential, and a shifted and truncated LJ potential have been examined. Although the bulk densities and surface tensions calculated for different model fluids are strongly affected by the LRC, the different potentials essentially lead to similar density values and similar surface tension values when the respective calculated values are compared on the basis of a reduced temperature scale.     

Vapor-liquid interfacial properties of rigid-linear Lennard-Jones chains

We have obtained the interfacial properties of short rigid-linear chains formed from tangentially bonded Lennard-Jones monomeric units from direct simulation of the vapour-liquid interface. The full long-range tails of the potential are accounted for by means of an improved version of the inhomogeneous long-range corrections of Janec?ek [J. Phys. Chem. B 110, 6264-6269 (2006)] proposed recently by MacDowell and Blas [J. Chem. Phys. 131, 074705 (2009)] valid for spherical as well as for rigid and flexible molecular systems. Three different model systems comprising of 3, 4, and 5 monomers per molecule are considered. The simulations are performed in the canonical ensemble, and the vapor-liquid interfacial tension is evaluated using the test-area method. In addition to the surface tension, we also obtain density profiles, coexistence densities, critical temperature and density, and interfacial thickness as functions of temperature, paying particular attention to the effect of the chain length and rigidity on these properties. According to our results, the main effect of increasing the chain length (at fixed temperature) is to sharpen the vapor-liquid interface and to increase the width of the biphasic coexistence region. As a result, the interfacial thickness decreases and the surface tension increases as the molecular chains get longer. The surface tension has been scaled by critical properties and represented as a function of the difference between coexistence densities relative to the critical density.     

Connection between the virial equation of state and physical clusters in a low density vapor

We carry out Monte Carlo simulations of physical Lennard-Jones and water clusters and show that the number of physical clusters in vapor is directly related to the virial equation of state. This relation holds at temperatures clearly below the critical temperatures, in other words, as long as the cluster-cluster interactions can be neglected--a typical assumption used in theories of nucleation. Above a certain threshold cluster size depending on temperature and interaction potential, the change in cluster work of formation can be calculated analytically with the recently proposed scaling law. The breakdown of the scaling law below the threshold sizes is accurately modeled with the low order virial coefficients. Our results indicate that high order virial coefficients can be analytically calculated from the lower order coefficients when the scaling law for cluster work of formation is valid. The scaling law also allows the calculation of the surface tension and equilibrium vapor density with computationally efficient simulations of physical clusters. Our calculated values are in good agreement with those obtained with other methods. We also present our results for the curvature dependent surface tension of water clusters.     

Long-range Lennard-Jones and electrostatic interactions in interfaces: application of the isotropic periodic sum method

Molecular dynamics (MD) simulations of heptane/vapor, hexadecane/vapor, water/vapor, hexadecane/water, and dipalmitoylphosphatidylcholine (DPPC) bilayers and monolayers are analyzed to determine the accuracy of treating long-range interactions in interfaces with the isotropic periodic sum (IPS) method. The method and cutoff (rc) dependences of surface tensions, density profiles, water dipole orientation, and electrostatic potential profiles are used as metrics. The water/vapor, heptane/vapor, and hexadecane/vapor interfaces are accurately and efficiently calculated with 2D IPS (rc=10 A). It is demonstrated that 3D IPS is not practical for any of the interfacial systems studied. However, the hybrid method PME/IPS [Particle Mesh Ewald for electrostatics and 3D IPS for Lennard-Jones (LJ) interactions] provides an efficient way to include both types of long-range forces in simulations of large liquid/vacuum and all liquid/liquid interfaces, including lipid monolayers and bilayers. A previously published pressure-based long-range LJ correction yields results similar to those of PME/IPS for liquid/liquid interfaces. The contributions to surface tension of LJ terms arising from interactions beyond 10 A range from 13 dyn/cm for the hexadecane/vapor interface to approximately 3 dyn/cm for hexadecane/water and DPPC bilayers and monolayers. Surface tensions of alkane/vapor, hexadecane/water, and DPPC monolayers based on the CHARMM lipid force fields agree very well with experiment, whereas surface tensions of the TIP3P and TIP4P-Ew water models underestimate experiment by 16 and 11 dyn/cm, respectively. Dipole potential drops (DeltaPsi) are less sensitive to long-range LJ interactions than surface tensions. However, DeltaPsi for the DPPC bilayer (845+/-3 mV proceeding from water to lipid) and water (547+/-2 mV for TIP4P-Ew and 521+/-3 mV for TIP3P) overestimate experiment by factors of 3 and 5, respectively, and represent expected deficiencies in nonpolarizable force fields.     

Effect of long-range interactions on the surface tension

The effect of long-range interactions on the surface tension at a liquid-gas interface was considered. An analytical expression for the correction to the surface tension for the cutoff of the particle interaction potential at the distance r _c was derived based on a step density profile. For the Lennard-Jones fluid, this correction was calculated numerically from the results of computer simulations of the density profiles. It was established that, in the vicinity of the triple point, the correction is as great as ～6% at the potential cutoff radius r _c=6.78 molecular diameters, a quantity insensitive to the form of the density profile in the interfacial layer.     

Surface tension of different sized single-component droplets,according to macroscopic data obtained using the lattice gas model and the critical droplet size during phase formation

Size dependences of the surface tension of spherical single-component droplets are calculated using equations of the lattice gas model for 19 compounds. Parameters of the model are found from experimental data on the surface tension of these compounds for a macroscopic planar surface. The chosen low-molecular compounds satisfy the law of corresponding states. To improve agreement with the experimental data, Lennard-Jones potential parameters are varied within 10% deviations. The surface tensions of different sized equilibrium droplets are calculated at elevated and lowered temperatures. It is found that the surface tension of droplets grows monotonically as the droplet size increases from zero to its bulk value. The droplet size R ₀ corresponding to zero surface tension corresponds to the critical size of the emergence of a new phase. The critical droplet sizes in the new phase of the considered compounds are estimated for the first time.     

Calculation of the surface tension of oxygen using molecular-dynamics simulations

The surface tension of oxygen at the liquid-vapor interface is calculated for the temperature range of 60-90 K using molecular-dynamics simulations and is shown to be within 1.0% error of experimental values for most of the temperatures studied. The potential used here is the same as in our previous study on liquid oxygen alone [S. D. Bembenek and B. M. Rice, J. Chem. Phys. 113, 2354 (2000)] and was optimized with an innovative statistical-mechanical method. The potential does not use a long-range cutoff nor a tail correction, which are usually considered necessary to obtain accurate values for the surface tension. We reason that the accuracy in surface tension is directly related to our parametrization method for the potential.     

Density, refractive index, interfacial tension, and viscosity of ionic liquids [EMIM][EtSO4], [EMIM][NTf2], [EMIM][N(CN)2], and [OMA][NTf2] in dependence on temperature at atmospheric pressure

The density, refractive index, interfacial tension, and viscosity of ionic liquids (ILs) [EMIM][EtSO 4] (1-ethyl-3-methylimidazolium ethylsulfate), [EMIM][NTf 2] (1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide), [EMIM][N(CN) 2] (1-ethyl-3-methylimidazolium dicyanimide), and [OMA][NTf 2] (trioctylmethylammonium bis(trifluoromethylsulfonyl)imide) were studied in dependence on temperature at atmospheric pressure both by conventional techniques and by surface light scattering (SLS). A vibrating tube densimeter was used for the measurement of density at temperatures from (273.15 to 363.15) K and the results have an expanded uncertainty ( k = 2) of +/-0.02%. Using an Abbe refractometer, the refractive index was measured for temperatures between (283.15 and 313.15) K with an expanded uncertainty ( k = 2) of about +/-0.0005. The interfacial tension was obtained from the pendant drop technique at a temperature of 293.15 K with an expanded uncertainty ( k = 2) of +/-1%. For higher and lower temperatures, the interfacial tension was estimated by an adequate prediction scheme based on the datum at 293.15 K and the temperature dependence of density. For the ILs studied within this work, at a first order approximation, the quantity directly accessible by the SLS technique was the ratio of surface tension to dynamic viscosity. By combining the experimental results of the SLS technique with density and interfacial tension from conventional techniques, the dynamic viscosity could be obtained for temperatures between (273.15 and 333.15) K with an estimated expanded uncertainty ( k = 2) of less than +/-3%. The measured density, refractive index, and viscosity are represented by interpolating expressions with differences between the experimental and calculated values that are comparable with but always smaller than the expanded uncertainties ( k = 2). Besides a comparison with the literature, the influence of structural variations on the thermophysical properties of the ILs is discussed in detail. The viscosities mostly agree with values reported in the literature within the combined estimated expanded uncertainties ( k = 2) of the measurements while our density and interfacial tension data differ by more than +/-1% and +/-5%.     

Parameterization and evaluation of a flexible water model

The thermodynamic, dielectric, and dynamic properties of a newly parameterized flexible water model are studied using molecular dynamics simulations. The potential function developed is based on the popular simple point charge (SPC) rigid model with the addition of appropriate harmonic and anharmonic energy terms for stretching and bending. Care was taken to account for the self-polarization and gas-phase monomer energy corrections during the parameterization, which have typically been ignored in past studies. The results indicate that an increased Lennard-Jones repulsive coefficient and slightly scaled partial charges are required when adding flexibility to the rigid model potential to reliably reproduce the experimental density, energy, and O ? O radial distribution function of water at 298 K and 1 atm. Analysis of the power spectrum derived from the H-velocity autocorrelation function allowed the water potential to be evaluated further and refined by adjusting the valence forces to fit the vibrational frequencies of the gas and liquid. Once a consistent set of parameters was determined, the static dielectric properties of the water model were calculated at two temperatures using the reaction field method to treat long-range forces and correlations. The dielectric constant of 75 ± 7 calculated at 300 K is in good agreement with the experimental value of 78.5. The Kirkwood g factor was also examined for temperature dependence and showed the correct increasing behavior with decreasing T. As a final check of the water potential, the free energies of solvation of a flexible water molecule and neon were predicted using thermodynamic perturbation methods. The calculated solvation energies of ?7.0 ± 0.8 for water and 2.7 ± 0.7 for neon are both consistent with the experimental values of ?6.3 and 2.7 kcal/mol. Comparisons are made throughout the study with the results of previous rigid and flexible model simulations. © 1995 by John Wiley & Sons, Inc.     

Surface tension of quantum fluids from molecular simulations

We present the first molecular simulations of the vapor-liquid surface tension of quantum liquids. The path integral formalism of Feynman was used to account for the quantum mechanical behavior of both the liquid and the vapor. A replica-data parallel algorithm was implemented to achieve good parallel performance of the simulation code on at least 32 processors. We have computed the surface tension and the vapor-liquid phase diagram of pure hydrogen over the temperature range 18-30 K and pure deuterium from 19 to 34 K. The simulation results for surface tension and vapor-liquid orthobaric densities are in very good agreement with experimental data. We have computed the interfacial properties of hydrogen-deuterium mixtures over the entire concentration range at 20.4 and 24 K. The calculated equilibrium compositions of the mixtures are in excellent agreement with experimental data. The computed mixture surface tension shows negative deviations from ideal solution behavior, in agreement with experimental data and predictions from Prigogine's theory. The magnitude of the deviations at 20.4 K are substantially larger from simulations and from theory than from experiments. We conclude that the experimentally measured mixture surface tension values are systematically too high. Analysis of the concentration profiles in the interfacial region shows that the nonideal behavior can be described entirely by segregation of H(2) to the interface, indicating that H(2) acts as a surfactant in H(2)-D(2) mixtures.     

Effect of flexibility on surface tension and coexisting densities of water

Molecular dynamics simulations of pure water at the liquid-vapor interface are performed using direct simulation of interfaces in a liquid slab geometry. The effect of intramolecular flexibility on coexisting densities and surface tension is analyzed. The dipole moment profile across the liquid-vapor interface shows different values for the liquid and vapor phases. The flexible model is a polarizable model. This effect is minor for liquid densities and is large for surface tension. The liquid densities increase from 2% at 300 K to 9% at 550 K when the force field is changed from a fully rigid simple point charge extended (SPCE) model to that of a fully flexible model with the same intermolecular interaction parameters. The increases in surface tension at both temperatures are around 11% and 36%, respectively. The calculated properties of the flexible models are closer to the experimental data than those of the rigid SPCE. The effect of the maximum number of reciprocal vectors (h(z) (max)) and the surface area on the calculated properties at 300 K is also analyzed. The coexiting densities are not sensitive to those variables. The surface tension fluctuates with h(z) (max) with an amplitude larger than 10 mN m(-1). The effect of using small interfacial areas is slightly larger than the error in the simulations.     

Computation of surface tensions using expanded ensemble simulations

A method for the direct simulation of the surface tension is examined. The technique is based on the thermodynamic route to the interfacial tension and makes use of the expanded ensemble simulation method for the calculation of the free energy difference between two inhomogeneous systems with the same number of particles, temperature, and volume, but different interfacial area. The method is completely general and suitable for systems with either continuous or discontinuous interactions. The adequacy of the expanded ensemble method is assessed by computing the interfacial tension of the planar vapor-liquid interface of Lennard-Jones, Lennard-Jones dimers, Gay-Berne, and square-well model fluids; in the latter, the interactions are discontinuous and the present method does not exhibit the asymmetry of other related methods, such as the test area. The expanded ensemble simulation results are compared with simulation data obtained from other techniques (mechanical and test area) with overall good agreement.     

Calculation of surface tension via area sampling

We examine the performance of several molecular simulation techniques aimed at evaluation of the surface tension through its thermodynamic definition. For all methods explored, the surface tension is calculated by approximating the change in Helmholtz free energy associated with a change in interfacial area through simulation of a liquid slab at constant particle number, volume, and temperature. The methods explored fall within three general classes: free-energy perturbation, the Bennett acceptance-ratio scheme, and the expanded ensemble technique. Calculations are performed for both the truncated Lennard-Jones and square-well fluids at select temperatures spaced along their respective liquid-vapor saturation lines. Overall, we find that Bennett and expanded ensemble approaches provide the best combination of accuracy and precision. All of the methods, when applied using sufficiently small area perturbation, generate equivalent results for the Lennard-Jones fluid. However, single-stage free-energy-perturbation methods and the closely related test-area technique recently introduced by Gloor et al. [J. Chem. Phys. 123, 134703 (2005)] generate surface tension values for the square-well fluid that are not consistent with those obtained from the more robust expanded ensemble and Bennett approaches, regardless of the size of the area perturbation. Single-stage perturbation methods fail also for the Lennard-Jones system when applied using large area perturbations. Here an analysis of phase-space overlap produces a quantitative explanation of the observed inaccuracy and shows that the satisfactory results obtained in these cases from the test-area method arise from a cancellation of errors that cannot be expected in general. We also briefly analyze the variation in method performance with respect to the adjustable parameters inherent to the techniques.     

Mean-field kinetic nucleation theory

A new semiphenomenological model of homogeneous vapor-liquid nucleation is proposed in which the cluster kinetics follows the "kinetic approach to nucleation" and the thermodynamic part is based on the revised Fisher droplet model with the mean-field argument for the cluster configuration integral. The theory is nonperturbative in a cluster size and as such is valid for all clusters down to monomers. It contains two surface tensions: macroscopic (planar) and microscopic. The latter is a temperature dependent quantity related to the vapor compressibility factor at saturation. For Lennard-Jones fluids the microscopic surface tension possesses a universal behavior with the parameters found from the mean-field density functional calculations. The theory is verified against nucleation experiments for argon, nitrogen, water, and mercury, demonstrating very good agreement with experimental data. Classical nucleation theory fails to predict experimental results when a critical cluster becomes small.     

On the applicability of Young–Laplace equation for nanoscale liquid drops

Debates continue on the applicability of the Young–Laplace equation for droplets, vapor bubbles and gas bubbles in nanoscale. It is more meaningful to find the error range of the Young–Laplace equation in nanoscale instead of making the judgement of its applicability. To do this, for seven liquid argon drops (containing 800, 1000, 1200, 1400, 1600, 1800, or 2000 particles, respectively) at T = 78 K we determined the radius of surface of tension R_s and the corresponding surface tension γ_s by molecular dynamics simulation based on the expressions of R_s and γ_s in terms of the pressure distribution for droplets. Compared with the two-phase pressure difference directly obtained by MD simulation, the results show that the absolute values of relative error of two-phase pressure difference given by the Young–Laplace equation are between 0.0008 and 0.027, and the surface tension of the argon droplet increases with increasing radius of surface of tension, which supports that the Tolman length of Lennard-Jones droplets is positive and that Lennard-Jones vapor bubbles is negative. Besides, the logic error in the deduction of the expressions of the radius and the surface tension of surface of tension, and in terms of the pressure distribution for liquid drops in a certain literature is corrected.     

A new method for computing surface tension using a drop of liquid

In the simulation of a liquid drop it is expensive to calculate the excess pressure and obtain the surface tension by the Laplace formula. We use the Kelvin formula which only requires the vapour density, or at most the virial pressure. Some results are given for a Lennard-Jones 12-6 fluid.