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.     

Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials

A novel test-area (TA) technique for the direct simulation of the interfacial tension of systems interacting through arbitrary intermolecular potentials is presented in this paper. The most commonly used method invokes the mechanical relation for the interfacial tension in terms of the tangential and normal components of the pressure tensor relative to the interface (the relation of Kirkwood and Buff [J. Chem. Phys. 17, 338 (1949)]). For particles interacting through discontinuous intermolecular potentials (e.g., hard-core fluids) this involves the determination of delta functions which are impractical to evaluate, particularly in the case of nonspherical molecules. By contrast we employ a thermodynamic route to determine the surface tension from a free-energy perturbation due to a test change in the surface area. There are important distinctions between our test-area approach and the computation of a free-energy difference of two (or more) systems with different interfacial areas (the method of Bennett [J. Comput. Phys. 22, 245 (1976)]), which can also be used to determine the surface tension. In order to demonstrate the adequacy of the method, the surface tension computed from test-area Monte Carlo (TAMC) simulations are compared with the data obtained with other techniques (e.g., mechanical and free-energy differences) for the vapor-liquid interface of Lennard-Jones and square-well fluids; the latter corresponds to a discontinuous potential which is difficult to treat with standard methods. Our thermodynamic test-area approach offers advantages over existing techniques of computational efficiency, ease of implementation, and generality. The TA method can easily be implemented within either Monte Carlo (TAMC) or molecular-dynamics (TAMD) algorithms for different types of interfaces (vapor-liquid, liquid-liquid, fluid-solid, etc.) of pure systems and mixtures consisting of complex polyatomic molecules.     

On interfacial tension calculation from the test-area methodology in the grand canonical ensemble

We propose the extension of the test-area methodology, originally proposed to evaluate the surface tension of planar fluid-fluid interfaces along a computer simulation in the canonical ensemble, to deal with the solid-fluid interfacial tension of systems adsorbed on slitlike pores using the grand canonical ensemble. In order to check the adequacy of the proposed extension, we apply the method for determining the density profiles and interfacial tension of spherical molecules adsorbed in slitlike pore with different pore sizes and solid-fluid dispersive energy parameters along the same simulation. We also calculate the solid-fluid interfacial tension using the original test-area method in the canonical ensemble. Agreement between the results obtained from both methods indicate that both methods are fully equivalent. The advantage of the new methodology is that allows to calculate simultaneously the density profiles and the amount of molecules adsorbed onto a slitlike pore, as well as the solid-fluid interfacial tension. This ensures that the chemical potential at which all properties are evaluated during the simulation is exactly the same since simulations can be performed in the grand canonical ensemble, mimicking the conditions at which the adsorption experiments are most usually carried out in the laboratory.     

Monte Carlo simulation methods for computing the wetting and drying properties of model systems

We introduce general Monte Carlo simulation methods for determining the wetting and drying properties of model systems. We employ an interface-potential-based approach in which the interfacial properties of a system are related to the surface excess free energy of a thin fluid film in contact with a surface. Two versions of this approach are explored: a "spreading" method focused on the growth of a thin liquid film from a surface in a mother vapor and a "drying" method focused on the growth of a thin vapor film from a surface in a mother liquid. The former provides a direct measure of the spreading coefficient while the latter provides an analogous drying coefficient. When coupled with an independent measure of the liquid-vapor surface tension, these coefficients enable one to compute the contact angle. We also show how one can combine information gathered from application of the spreading and drying methods at a common state point to obtain direct measures of the contact angle and liquid-vapor surface tension. The computational strategies introduced here are applied to two model systems. One includes a monatomic Lennard-Jones fluid that interacts with a structureless substrate via a long-ranged substrate potential. The second model contains a monatomic Lennard-Jones fluid that interacts with an atomistically detailed substrate via a short-ranged potential. Expanded ensemble techniques are coupled with the interface potential approach to compile the temperature- and substrate strength-dependence of various interfacial properties for these systems. Overall, we find that the approach pursued here provides an efficient and precise means to calculate the wetting and drying properties of model systems.     

Universal scaling behaviour of surface tension of molecular chains

We use and extend the universal relationship recently proposed by Galliero [G. Galliero, J. Chem. Phys. 133, 074705 (2010)], based on a combination of the corresponding-states principle of Guggenheim [E. A. Guggenheim, J. Chem. Phys. 13, 253 (1945)] and the parachor approach of Macleod [J. Macleod, Trans. Faraday Soc. 19, 38 (1923)], to predict the vapour-liquid surface tension of fully flexible chainlike Lennard-Jones molecules. In the original study of Galliero, the reduced surface tension of short-chain molecules formed by up to five monomers is expressed as a unique function of the difference between the liquid and vapour coexistence densities. In this work, we extend the applicability of the recipe and demonstrate that it is also valid for predicting the surface tension of two different chainlike molecular models, namely, linear tangent chains that interact through the Lennard-Jones intermolecular potential and fully flexible chains formed by spherical segments interacting through the square-well potential. Computer simulation data for vapour-liquid surface tension of fully flexible and rigid linear Lennard-Jones, and fluid flexible square-well chains is taken from our previous works. Our results indicate that the universal scaling relationship is able to correlate short- and long-chain molecules with different degrees of flexibility and interacting through different intermolecular potentials.     

Calculation of the surface tension and pressure components from a non-exponential perturbation method of the thermodynamic route

Surface tension is probably the most important interfacial property and a large number of techniques have been devoted to its calculation. Usually, this calculation is carried out using mechanical or thermodynamic definitions. The mechanical route uses an arbitrary choice to affect the contribution of the pairwise force. To overcome this arbitrariness, a thermodynamic route based on the area perturbation (test-area (TA) method) has been developed for the calculation of surface tension. The volume perturbation (VP) method provides an original route to compute the components of the pressure tensor. These two routes are developed from the perturbation theory leading to working expressions using exponential averages of energy. The use of exponential averages makes the calculation strongly dependent on the occurrence of low values of ΔU. Additionally, the decomposition of the energy to obtain local surface tension is nontrivial. From the explicit derivation of the partition function the exponential average is avoided providing an interesting alternative to TA, VP, and mechanical methods. To make a consistent comparison, we study the profiles of the surface tension along the direction normal to the surface for the different definitions and techniques in the cases of liquid-vapor interfaces of acids gases, binary, and apolar systems.     

A unified methodological framework for the simulation of nonisothermal ensembles

A general framework is developed for the simulation of nonisothermal statistical-mechanical ensembles. This framework is intended to synthesize the formulation of advanced Monte Carlo simulation methods such as multihistogram reweighting, replica-exchange methods, and expanded ensemble techniques so that they can be applied to different nonisothermal ensembles. Using Lennard-Jones systems as test cases, novel implementations of these methods are demonstrated with different ensembles including the microcanonical, isobaric-isoenthalpic, and isobaric-semigrand ensembles. In particular, it is shown that the use of multiensemble methods allows the efficient simulation of microcanonical density of states, entropies, vapor-liquid and solid-liquid equilibrium for pure component systems, and fluid-phase coexistence for binary mixtures. In these applications, comparisons are also presented that highlight the advantages of the proposed multiensemble implementations over alternative methods used before.     

Calculating thermodynamic properties from perturbation theory: II. An analytic representation for the square-well chain fluid

An analytic representation of thermodynamic properties of the freely jointed square-well chain fluid is developed based on the thermodynamic perturbation theory of Barker–Henderson, Zhang and Weitheim. By using a real function expression for the radial distribution function and incorporating structural information for square-well monomer of TPT1 model, an analytic expression for the Helmholtz energy of square-well chain fluid is expanded from Zhang’s analytic expressions for thermodynamic properties of square-well monomer. The expression leads to good predictions of the compressibility factor, residual internal energy and constant-volume heat capacity for 4-mer, 8-mer and 16-mer square-well fluids when compared with the Monte Carlo (MC) simulation results. The incorporating structural information for square-well dimer of TPT-D model is also calculated. To obtain the constant-volume heat capacity needed, NVT MC simulations were performed.     

Simulation of chemical potentials and phase equilibria in two- and three-dimensional square-well fluids: finite size effects

We study the simulation cell size dependence of chemical potential isotherms in subcritical square-well fluids by means of series of canonical ensemble Monte Carlo simulations with increasing numbers of particles, for both three-dimensional bulk systems and two-dimensional planar layers, using Widom-like particle insertion methods. By estimating the corresponding vapor/liquid coexistence densities using a Maxwell-like equal area rule for the subcritical chemical potential isotherms, we are able to study the influence of system size not only on chemical potentials but also on the coexistence properties. The chemical potential versus density isotherms show van der Waals-like loops in the subcritical vapor/liquid coexistence range that exhibit distinct finite size effects for both two- and three-dimensional fluids. Generally, in agreement with recent findings for related studies of Lennard-Jones fluids, the loops shrink with increasing number of particles. In contrast to the subcritical isotherms themselves, the equilibrium vapor/liquid densities show only a weak system size dependence and agree quantitatively with the best-known literature values for three-dimensional fluids. This allows our approach to be used to accurately predict the phase coexistence properties. Our resulting phase equilibrium results for two-dimensional square-well fluids are new. Knowledge concerning finite size effects of square-well systems is important not only for the simulation of thermodynamic properties of simple fluids, but also for the simulation of models of more complex fluids (such as aqueous or polymer fluids) involving square-well interactions.     

Oscillatory surface tension due to finite-size effects

The simulation results of surface tension at the liquid-vapor interface are presented for fluids interacting with Lennard Jones and square-well potentials. From the simulation of liquids we have reported [M. González-Melchor et al., J. Chem. Phys. 122, 4503 (2005)] that the components of pressure tensor in parallelepiped boxes are not the same when periodic boundary conditions and small transversal areas are used. This fact creates an artificial oscillatory stress anisotropy in the system with even negative values. By doing direct simulations of interfaces we show in this work that surface tension has also an oscillatory decay at small surface areas; this behavior is opposite to the monotonic decay reported previously for the Lennard Jones fluid. It is shown that for small surface areas, the surface tension of the square-well potential artificially takes negative values and even increases with temperature. The calculated surface tension using a direct simulation of interfaces might have two contributions: one from finite-size effects of interfacial areas due to box geometry and another from the interface. Thus, it is difficult to evaluate the true surface tension of an interface when small surface areas are used. Care has to be taken to use the direct simulation method of interfaces to evaluate the predicted surface tension as a function of interfacial area from capillary-wave theory. The oscillations of surface tension decay faster at temperatures close to the critical point. It is also discussed that a surface area does not show any important effect on coexisting densities, making this method reliable to calculate bulk coexisting properties using small systems.     

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.     

Second order of perturbation theory correction to the radial distribution function of a liquid with the interaction potential of hard spheres plus a square-well

A liquid with the interaction potential of hard spheres plus a square-well is analyzed using the Monte-Carlo technique. Numerical results for the perturbation theory series over a square-well potential are obtained in the form of the Barker and Henderson discrete representation. Approximating expressions for the correction to a liquid radial distribution function in the second order of perturbation theory are presented. The obtained results allow us to define this correction with a root-mean-square deviation of about 0.007. It is shown that the given approach provides a complete calculation in the second order of perturbation theory, and also the determination of the third order correction to the free energy for a liquid interacting with the potential of the Lennard-Jones type.     

Surface tension of a Lennard-Jones liquid under supersaturation

A formally exact Kirkwood-Buff virial formula for the surface tension of a supersaturated interface is derived. A modified Gibbs ensemble method is given that allows the creation of interacting supersaturated phases of equal chemical potential, and which enables the Kirkwood-Buff formula to be applied. The methods are tested by Monte Carlo simulation of a supersaturated Lennard-Jones fluid with a planar liquid-vapor interface. The Kirkwood-Buff results for the supersaturated surface tension are found to be in reasonable agreement with new results obtained here using the recently developed, formally exact, ghost interface method, [M. P. Moody and P. Attard, J. Chem. Phys., 2002, 117, 6705]. The surface tension is obtained as a function of supersaturation at four temperatures, and it is found to decrease with increasing supersaturation, and to vanish at the vapor spinodal. The relevance of the present results to the nucleation of droplets in a supersaturated vapor is discussed.     

Effective potential approach to bulk thermodynamic properties and surface tension of molecular fluids I. Single component n-alkane systems

The aim of this work is to develop spherically symmetric effective potentials allowing bulk thermodynamic properties and surface tension of molecular fluids to be predicted semiempirically by the use of statistical mechanical methods. Application is made to the straight chain alkane fluids from methane to decane. An effective Lennard-Jones potential is generated with temperature-dependent parameters fitted to the critical temperature and pressure and to Pitzer's acentric factor. Insertion of this potential into the generalised van der Waals (GvdW) density functional theory yields bulk properties in good agreement with experiments. The surface tension is overestimated for the longer alkane chains. In order to account for the surface tension, an independently adjustable attractive range of interaction is required and obtained through the use of square-well potentials chosen so as to leave the bulk thermodynamics unaltered while the attractive range is fitted to the surface tension at a single temperature. The GvdW theory, which includes binding energy, entropic and profile shape contributions, then generates surface tension estimates that are of good accuracy over the full range of available experimental data. It appears that, given a sufficiently flexible form, effective potentials combined with simple statistical mechanical theory can reproduce both bulk and non-uniform fluid data of great variety in an insighful and practically useful way.     

Isothermalisobaric molecular dynamics simulation of diatomic liquids and their mixtures

Isothermalisobaric ensemble molecular dynamics simulations have been performed for systems of two-center Lennard-Jones molecules for pure fluids as well as binary mixtures. The results obtained from various ensemble averages have been compared for pure fluid simulations of Lennard-Jones model diatomic fluids. The excess enthalpy and excess volume of three equimolar mixtures (argonnitrogen, argonoxygen, and nitrogenoxygen) have been calculated and compared with values obtained from previous NVT molecular dynamics and perturbation theory. Pair distribution functions obtained from various methods are compared for the equimolar mixture of nitrogen and oxygen and used to study the effect of attractive forces on the local structure. For four other systems (argonethane, methaneethane, carbon disulfidecarbon tetrachloride, and carbon disulfidetetrachloroethylene), excess enthalpies and excess volumes from NPT simulations are used to test the ability of perturbation theory to predict these properties and are also compared with experimental data. The comparison of simulation and experiment clearly shows the need to improve the available parameters for cross interactions in binary mixtures.     

Surface tension and vapor-liquid phase coexistence of confined square-well fluid

Phase equilibria of a square-well fluid in planar slit pores with varying slit width are investigated by applying the grand-canonical transition-matrix Monte Carlo (GC-TMMC) with the histogram-reweighting method. The wall-fluid interaction strength was varied from repulsive to attractive such that it is greater than the fluid-fluid interaction strength. The nature of the phase coexistence envelope is in agreement with that given in literature. The surface tension of the vapor-liquid interface is calculated via molecular dynamics simulations. GC-TMMC with finite size scaling is also used to calculate the surface tension. The results from molecular dynamics and GC-TMMC methods are in very good mutual agreement. The vapor-liquid surface tension, under confinement, was found to be lower than the bulk surface tension. However, with the increase of the slit width the surface tension increases. For the case of a square-well fluid in an attractive planar slit pore, the vapor-liquid surface tension exhibits a maximum with respect to wall-fluid interaction energy. We also report estimates of critical properties of confined fluids via the rectilinear diameter approach.     

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.     

Thin-thick surface phase coexistence and boundary tension of the square-well fluid on a weak attractive surface

Prewetting transition is studied for the square-well fluid of attractive-well diameter lambda(ff)sigma(ff)=1.5 in the presence of a homogeneous surface modeled by the square-well potential of attractive well from 0.8sigma(ff) to 1.8sigma(ff). We investigate surface phase coexistence of thin-thick film transition using grand-canonical transition matrix Monte Carlo (GC-TMMC) and histogram reweighting techniques. Molecular dynamics (MD) and GC-TMMC are utilized to predict the properties of the fluid for various surface fluid affinities. Occurrences of prewetting transition with the variation of surface affinity are observed for a domain of reduced temperature from T(*)=0.62 to 0.75. We have used MD and GC-TMMC+finite size scaling (FSS) simulations to calculate the boundary tension as a function of temperature as well as surface affinity. Boundary tensions via MD and GC-TMMC+FSS methods are in good agreement. The boundary tension increases with the decrease of wall-fluid affinity. Prewetting critical properties are calculated using rectilinear diameter approach and scaling analysis. We found that critical temperature and density increase with the decrease of wall-fluid affinity.     

Density functional theory of inhomogeneous liquids. I. The liquid-vapor interface in Lennard-Jones fluids

A simple model is proposed for the direct correlation function (DCF) for simple fluids consisting of a hard-core contribution, a simple parametrized core correction, and a mean-field tail. The model requires as input only the free energy of the homogeneous fluid, obtained, e.g., from thermodynamic perturbation theory. Comparison to the DCF obtained from simulation of a Lennard-Jones fluid shows this to be a surprisingly good approximation for a wide range of densities. The model is used to construct a density functional theory for inhomogeneous fluids which is applied to the problem of calculating the surface tension of the liquid-vapor interface. The numerical values found are in good agreement with simulation.     

Vapor-liquid equilibrium and critical behavior of the square-well fluid of variable range: a theoretical study

The vapor-liquid phase behavior and the critical behavior of the square-well (SW) fluid are investigated as a function of the interaction range, lambdain [1.25, 3], by means of the self-consistent Ornstein-Zernike approximation (SCOZA) and analytical equations of state based on a perturbation theory [A. L. Benavides and F. del Rio, Mol. Phys. 68, 983 (1989); A. Gil-Villegas, F. del Rio, and A. L. Benavides, Fluid Phase Equilib. 119, 97 (1996)]. For this purpose the SCOZA, which has been restricted up to now to a few model systems, has been generalized to hard-core systems with arbitrary interaction potentials requiring a fully numerical solution of an integro-partial differential equation. Both approaches, in general, describe well the liquid-vapor phase diagram of the square-well fluid when compared with simulation data. SCOZA yields very precise predictions for the coexistence curves in the case of long ranged SW interaction (lambda>1.5), and the perturbation theory is able to predict the binodal curves and the saturated pressures, for all interaction ranges considered if one stays away from the critical region. In all cases, the SCOZA gives very good predictions for the critical temperatures and the critical pressures, while the perturbation theory approach tends to slightly overestimate these quantities. Furthermore, we propose analytical expressions for the critical temperatures and pressures as a function of the square-well range.