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.     

Determination of the melting point of hard spheres from direct coexistence simulation methods

We consider the computation of the coexistence pressure of the liquid-solid transition of a system of hard spheres from direct simulation of the inhomogeneous system formed from liquid and solid phases separated by an interface. Monte Carlo simulations of the interfacial system are performed in three different ensembles. In a first approach, a series of simulations is carried out in the isothermal-isobaric ensemble, where the solid is allowed to relax to its equilibrium crystalline structure, thus avoiding the appearance of artificial stress in the system. Here, the total volume of the system fluctuates due to changes in the three dimensions of the simulation box. In a second approach, we consider simulations of the inhomogeneous system in an isothermal-isobaric ensemble where the normal pressure, as well as the area of the (planar) fluid-solid interface, are kept constant. Now, the total volume of the system fluctuates due to changes in the longitudinal dimension of the simulation box. In both approaches, the coexistence pressure is estimated by monitoring the evolution of the density along several simulations carried out at different pressures. Both routes are seen to provide consistent values of the fluid-solid coexistence pressure, p=11.54(4)k(B)T/sigma(3), which indicates that the error introduced by the use of the standard constant-pressure ensemble for this particular problem is small, provided the systems are sufficiently large. An additional simulation of the interfacial system is conducted in a canonical ensemble where the dimensions of the simulation box are allowed to change subject to the constraint that the total volume is kept fixed. In this approach, the coexistence pressure corresponds to the normal component of the pressure tensor, which can be computed as an appropriate ensemble average in a single simulation. This route yields a value of p=11.54(4)k(B)T/sigma(3). We conclude that the results obtained for the coexistence pressure from direct simulations of the liquid and solid phases in coexistence using different ensembles are mutually consistent and are in excellent agreement with the values obtained from free energy calculations.     

An exact formalism to study the thermodynamic properties of hard-sphere systems under spherical confinement

This paper presents a modified grand canonical ensemble which provides a new simple and efficient scheme to study few-body fluid-like inhomogeneous systems under confinement. The new formalism is implemented to investigate the exact thermodynamic properties of a hard sphere (HS) fluid-like system with up to three particles confined in a spherical cavity. In addition, the partition function of this system was used to analyze the surface thermodynamic properties of the many-HS system and to derive the exact curvature dependence of both the surface tension and adsorption in powers of the density. The expressions for the surface tension and the adsorption were also obtained for the many-HS system outside of a fixed hard spherical object. We used these results to derive the dependence of the fluid-substrate Tolman length up to first order in density.     

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.     

Stress anisotropy induced by periodic boundary conditions

Finite size effects due to periodic boundary conditions are investigated using computer simulations in the canonical ensemble. We study liquids with densities corresponding to typical liquid coexistence densities, and temperatures between the triple and critical points. The components of the pressure tensor are computed in order to analyze the finite size effects arising from the size and geometry of the simulation box. Two different box geometries are considered: cubic and parallelepiped. As expected the pressure tensor is isotropic in cubic boxes, but it becomes anisotropic for small noncubic boxes. We argue this is the origin of the anomalous behavior observed recently in the computation of the surface tension of liquid-vapor interfaces. Otherwise, we find that the bulk pressure is sensitive to the box geometry when small simulation boxes are considered. These observations are general and independent of the model liquid considered. We report results for liquids interacting through short range forces, square well and Lennard-Jones, and also long range Coulombic interactions. The effect that small surface areas have on the surface tension is discussed, and some preliminary results at the liquid vapor-interface for the square well potential are given.     

Finite-size effects in dissipative particle dynamics simulations

We have performed dissipative particle dynamics (DPD) simulations to evaluate the effect that finite size of transversal area has on stress anisotropy and interfacial tension. The simulations were carried out in one phase and two phases in parallelepiped cells. In one-phase simulations there is no finite-size effect on stress anisotropy when the simulation is performed using repulsive forces. However, an oscillatory function of stress anisotropy is found for attractive-repulsive interactions. In the case of liquid-liquid interfaces with repulsive interaction between molecules, there is only a small effect of surface area on interfacial tension when the simulations are performed using the Monte Carlo method at constant temperature and normal pressure. An important but artificial finite-size effect of interfacial area on surface tension is found in simulations in the canonical ensemble. Reliable results of interfacial tension from DPD simulations can be obtained using small systems, less than 2000 particles, when they interact exclusively with repulsive forces.     

Partial rigid-body dynamics in NPT, NPAT and NPgammaT ensembles for proteins and membranes

A partial rigid-body method of molecular dynamics simulations for proteins and membranes is presented. In this method, the symplectic integrator for rigid bodies is combined with the equations of motion for the NPT ensemble. The standard NPT ensemble is extended to the membrane-specific ensembles, the NPAT (constant normal pressure and lateral surface area of membranes and constant temperature) and NPgammaT (constant normal pressure and lateral surface tension of membranes and constant temperature) ensembles. By more than 30-ns simulations of aqueous proteins and hydrated lipid bilayers, the results of the partial rigid-body method demonstrated excellent conservation of total energy and consistent behavior with the traditional constraint method in terms of structural distribution and fluctuation of proteins and lipids. The efficient implementation of the partial rigid-body method in parallel computation is presented, which is shown to work well in large-scale molecular dynamics simulations.     

Multiple histogram reweighting method for the surface tension calculation

The multiple histogram reweighting method takes advantage of calculating ensemble averages over a range of thermodynamic conditions without performing a molecular simulation at each thermodynamic point. We show that this method can easily be extended to the calculation of the surface tension. We develop a new methodology called multiple histogram reweighting with slab decomposition based on the decomposition of the system into slabs along the direction normal to the interface. The surface tension is then calculated from local values of the chemical potential and of the configurational energy using Monte Carlo (MC) simulations. We show that this methodology gives surface tension values in excellent agreement with experiments and with standard NVT MC simulations in the case of the liquid-vapor interface of carbon dioxide.     

The nature of the calculation of the pressure in molecular simulations of continuous models from volume perturbations

We consider some fundamental aspects of the calculation of the pressure from simulations by performing volume perturbations. The method, initially proposed for hard-core potentials by Eppenga and Frenkel [Mol. Phys.52, 1303 (1984)] and then extended to continuous potentials by Harismiadis et al. [J. Chem. Phys. 105, 8469 (1996)], is based on the numerical estimate of the change in Helmholtz free energy associated with the perturbation which, in turn, can be expressed as an ensemble average of the corresponding Boltzmann factor. The approach can be easily generalized to the calculation of components of the pressure tensor and also to ensembles other than the canonical ensemble. The accuracy of the method is assessed by comparing simulation results obtained from the volume-perturbation route with those obtained from the usual virial expression for several prototype fluid models. Monte Carlo simulation data are reported for bulk fluids and for inhomogeneous systems containing a vapor-liquid interface.     

The interfacial tension and phase diagram of the Widom-Rowlinson mixture via Monte Carlo simulations

Results of Monte Carlo simulations are reported for the interfacial tension between two fluid phases in a binary mixture of penetrable spheres in which molecular pairs of like species do not interact, while those of unlike species interact as hard spheres. Semigrand canonical ensemble Monte Carlo simulations in a cubic cell with periodic boundary conditions are used to obtain histograms for various system sizes at various densities. At a given density, the interfacial tension and compositions of coexisting phases for an infinite system are evaluated via histogram analysis combined with finite-size scaling. The density dependence of the interfacial tension and phase diagram for an infinite system are thus obtained. The simulated behavior of the interfacial tension close to the critical density corroborates previous suggestions that the model belongs to the three-dimensional Ising universality class.     

Molecular simulations of confined liquids: an alternative to the grand canonical Monte Carlo simulations

Commonly, the confinement effects are studied from the grand canonical Monte Carlo (GCMC) simulations from the computation of the density of liquid in the confined phase. The GCMC modeling and chemical potential (μ) calculations are based on the insertion/deletion of the real and ghost particle, respectively. At high density, i.e., at high pressure or low temperature, the insertions fail from the Widom insertions while the performing methods as expanded method or perturbation approach are not efficient to treat the large and complex molecules. To overcome this problem we use a simple and efficient method to compute the liquid's density in the confined medium. This method does not require the precalculation of μ and is an alternative to the GCMC simulations. From the isothermal-isosurface-isobaric statistical ensemble we consider the explicit framework/liquid external interface to model an explicit liquid's reservoir. In this procedure only the liquid molecules undergo the volume changes while the volume of the framework is kept constant. Therefore, this method is described in the Np(n)AV(f)T statistical ensemble, where N is the number of particles, p(n) is the normal pressure, V(f) is the volume of framework, A is the surface of the solid/fluid interface, and T is the temperature. This approach is applied and validated from the computation of the density of the methanol and water confined in the mesoporous cylindrical silica nanopores and the MIL-53(Cr) metal organic framework type, respectively.     

Interfacial tension and interfacial profiles: an equation-of-state approach

A quasi-thermodynamic approach of inhomogeneous systems is used for modeling the fluid-fluid interface. It is based on the recently introduced QCHB (quasi-chemical hydrogen bonding) equation-of-state model of fluids and their mixtures, which is used for the estimation of the Helmholtz free energy density difference, Deltapsi(0), between the system with interface and another system of the same constitution but without interface. Consistent expressions for the interfacial tension and interfacial profiles for various properties are presented. The interfacial tension is proportional to the integral of Deltapsi(0) along the full height of the system, the proportionality constant being equal to 1, when no density gradient contributions are taken into consideration, 2, when the Cahn-Hilliard approximation is adopted, and 4, when the full density gradient contributions are taken into consideration. A satisfactory agreement is obtained between experimental and calculated surface tensions. Extension of the approach to mixtures is examined along with the associated problems for the numerical calculations of the interfacial profiles. A new equation is derived for the chemical potentials in the interfacial region, which facilitates very much the calculation of the composition profiles across the interface.     

Application of kinetic Monte Carlo method to equilibrium systems: vapour-liquid equilibria

Kinetic Monte Carlo (kMC) simulations were carried out to describe the vapour-liquid equilibria of argon at various temperatures. This paper aims to demonstrate the potential of the kMC technique in the analysis of equilibrium systems and its advantages over the traditional Monte Carlo method, which is based on the Metropolis algorithm. The key feature of the kMC is the absence of discarded trial moves of molecules, which ensures larger number of configurations that are collected for time averaging. Consequently, the kMC technique results in significantly fewer errors for the same number of Monte Carlo steps, especially when the fluid is rarefied. An additional advantage of the kMC is that the relative displacement probability of molecules is significantly larger in rarefied regions, which results in a more efficient sampling. This provides a more reliable determination of the vapour phase pressure and density in case of non-uniform density distributions, such as the vapour-liquid interface or a fluid adsorbed on an open surface. We performed kMC simulations in a canonical ensemble, with a liquid slab in the middle of the simulation box to model two vapour-liquid interfaces. A number of thermodynamic properties such as the pressure, density, heat of evaporation and the surface tension were reliably determined as time averages.     

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.     

Phase separation of a model binary polymer solution in an external field

The phase separation of a simple binary mixture of incompatible linear polymers in solution is investigated using an extension of the sedimentation equilibrium method, whereby the osmotic pressure of the mixture is extracted from the density profiles of the inhomogeneous mixture in a gravitational field. In Monte Carlo simulations the field can be tuned to induce significant inhomogeneity, while keeping the density profiles sufficiently smooth for the macroscopic condition of hydrostatic equilibrium to remain applicable. The method is applied here for a simplified model of ideal but mutually avoiding polymers, which readily phase separate at relatively low densities. The Monte Carlo data are interpreted with the help of an approximate bulk phase diagram calculated from a simple, second-order virial coefficient theory. By derivation of effective potentials between polymer centers of mass, the binary mixture of polymers is coarse-grained to a "soft colloid" picture reminiscent of the Widom-Rowlinson model for incompatible atomic mixtures. This approach significantly speeds up the simulations and accurately reproduces the behavior of the full monomer resolved model.     

Capillary waves at the liquid-vapor interface and the surface tension of water

Capillary waves occurring at the liquid-vapor interface of water are studied using molecular dynamics simulations. In addition, the surface tension, determined thermodynamically from the difference in the normal and tangential pressure at the liquid-vapor interface, is compared for a number of standard three- and four-point water models. We study four three-point models (SPC/E, TIP3P, TIP3P-CHARMM, and TIP3P-Ew) and two four-point models (TIP4P and TIP4P-Ew). All of the models examined underestimate the surface tension; the TIP4P-Ew model comes closest to reproducing the experimental data. The surface tension can also be determined from the amplitude of capillary waves at the liquid-vapor interface by varying the surface area of the interface. The surface tensions determined from the amplitude of the logarithmic divergence of the capillary interfacial width and from the traditional thermodynamic method agree only if the density profile is fitted to an error function instead of a hyperbolic tangent function.     

Molecular dynamics simulation of nanodroplet spreading enhanced by linear surfactants

We utilize molecular dynamics simulations to probe the surfactant-mediated spreading of a Lennard-Jones liquid droplet on a solid surface. The surfactants are linear hexamers that are insoluble in the liquid and reduce the surface tension of the liquid-vapor interface. We study how the interaction of the surfactant hexamers with the solid substrate influences spreading, as well as the dependence of spreading on surfactant concentration. We find that the spreading speed is strongly influenced by the attraction of the hydrophobic surfactant tail to the solid surface. When this attraction is sufficiently strong, surfactant molecules partition to the liquid-solid interface and facilitate spreading. This partitioning can lead to an inhomogeneous distribution of surfactant over the liquid-vapor interface, which could drive the Marangoni convection. We also observe that the surfactant molecules can assemble into micelles on the solid surface. The repulsion between micelles at the liquid-solid interface can lead to break-off and migration of the micelles from the liquid-solid to the gas-solid interface and spreading is facilitated in this way. Our model system contains features that are believed to underlie superspreading in experimental studies of droplet spreading.     

Enabling grand‐canonical Monte Carlo: Extending the flexibility of GROMACS through the GromPy python interface module

We report on a python interface to the GROMACS molecular simulation package, GromPy (available at https://github.com/GromPy ). This application programming interface (API) uses the ctypes python module that allows function calls to shared libraries, for example, written in C. To the best of our knowledge, this is the first reported interface to the GROMACS library that uses direct library calls. GromPy can be used for extending the current GROMACS simulation and analysis modes. In this work, we demonstrate that the interface enables hybrid Monte‐Carlo/molecular dynamics (MD) simulations in the grand‐canonical ensemble, a simulation mode that is currently not implemented in GROMACS. For this application, the interplay between GromPy and GROMACS requires only minor modifications of the GROMACS source code, not affecting the operation, efficiency, and performance of the GROMACS applications. We validate the grand‐canonical application against MD in the canonical ensemble by comparison of equations of state. The results of the grand‐canonical simulations are in complete agreement with MD in the canonical ensemble. The python overhead of the grand‐canonical scheme is only minimal. © 2012 Wiley Periodicals, Inc.     

Monte Carlo simulations of cylinder-forming ABC triblock terpolymer thin films

We systematically study the cylinder-forming ABC triblock terpolymer thin films using canonical ensemble Monte Carlo simulations. The simulated annealing procedure is applied to the self-assembling process. By judicious choice of the system dimensions, we elaborately investigate the effect of film thickness on the orientation of the cylinders. This confined triblock terpolymer system exhibits different phase behavior under the weak and strong surface fields. In addition, we also investigate the ensemble-averaged chain orientations and relative density profiles.     

A simulation method for the calculation of chemical potentials in small, inhomogeneous, and dense systems

We present a modification of the gauge cell Monte Carlo simulation method [A. V. Neimark and A. Vishnyakov, Phys. Rev. E 62, 4611 (2000)] designed for chemical potential calculations in small confined inhomogeneous systems. To measure the chemical potential, the system under study is set in chemical equilibrium with the gauge cell, which represents a finite volume reservoir of ideal particles. The system and the gauge cell are immersed into the thermal bath of a given temperature. The size of the gauge cell controls the level of density fluctuations in the system. The chemical potential is rigorously calculated from the equilibrium distribution of particles between the system cell and the gauge cell and does not depend on the gauge cell size. This scheme, which we call a mesoscopic canonical ensemble, bridges the gap between the canonical and the grand canonical ensembles, which are known to be inconsistent for small systems. The ideal gas gauge cell method is illustrated with Monte Carlo simulations of Lennard-Jones fluid confined to spherical pores of different sizes. Special attention is paid to the case of extreme confinement of several molecular diameters in cross section where the inconsistency between the canonical ensemble and the grand canonical ensemble is most pronounced. For sufficiently large systems, the chemical potential can be reliably determined from the mean density in the gauge cell as it was implied in the original gauge cell method. The method is applied to study the transition from supercritical adsorption to subcritical capillary condensation, which is observed in nanoporous materials as the pore size increases.