Efficiency Comparison of Single- and Multiple-Macrostate Grand Canonical Ensemble Transition-Matrix Monte Carlo Simulations.

Recent interest in parallelizing flat-histogram transition-matrix Monte Carlo simulations in the grand canonical ensemble, due to its demonstrated effectiveness in studying phase behavior, self-assembly and adsorption, has led to the most extreme case of single-macrostate simulations, where each macrostate is simulated independently with ghost particle insertions and deletions. Despite their use in several studies, no efficiency comparisons of these single-macrostate simulations have been made with multiple-macrostate simulations. We show that multiple-macrostate simulations are up to 3 orders of magnitude more efficient than single-macrostate simulations, which demonstrates the remarkable efficiency of flat-histogram biased insertions and deletions, even with low acceptance probabilities. Efficiency comparisons were made for supercritical fluids and vapor-liquid equilibrium of bulk Lennard-Jones and a three-site water model, self-assembling patchy trimer particles and adsorption of a Lennard-Jones fluid confined in a purely repulsive porous network, using the open source simulation toolkit FEASST. By directly comparing with a variety of Monte Carlo trial move sets, this efficiency loss in single-macrostate simulations is attributed to three related reasons. First, ghost particle insertions and deletions in single-macrostate simulations incur the same computational expense as grand canonical ensemble trials in multiple-macrostate simulations, yet ghost trials do not reap the sampling benefit from propagating the Markov chain to a new microstate. Second, single-macrostate simulations lack macrostate change trials that are biased by the self-consistently converging relative macrostate probability, which is a major component of flat histogram simulations. Third, limiting a Markov chain to a single macrostate reduces sampling possibilities. Existing parallelization methods for multiple-macrostate flat-histogram simulations are shown to be more efficient than parallel single-macrostate simulations by approximately an order of magnitude or more in all systems investigated.

Coupled Monte Carlo and Molecular Dynamics Simulations on Interfacial Properties of Antifouling Polymer Membranes.

We use molecular simulation to study the wetting behavior of antifouling polymer-tethered membranes. We obtain the interfacial properties (e.g., contact angle) of water at various temperatures for five polymer membranes, including a base polysulfone (PSF) membrane and four other PSF membranes grafted with antifouling polymers (two poly(ethylene glycol) (PEG) tethers and two zwitterionic tethers). We implement a coupled Monte Carlo (MC)/molecular dynamics (MD) approach to determine the interface potentials of water on the membrane surfaces in an efficient manner. Within this method, short MC and MD simulations are performed in cycles to collect the surface excess free energy of a thin water film on polymer membrane surfaces. Simulation results show that the grafting of zwitterionic tethers provides a more significant enhancement in the hydrophilicity of the PSF membrane than that of the PEG tethers. Water completely wets the surface of zwitterionic polymer membranes.

Improving the efficiency of Monte Carlo simulations of ions using expanded grand canonical ensembles.

While ionic liquids have promising applications as industrial solvents, predicting their fluid phase properties and coexistence remains a challenge. Grand canonical Monte Carlo simulation is an effective method for such predictions, but equilibration is hampered by the apparent requirement to insert and delete neutral sets of ions simultaneously in order to maintain charge neutrality. For relatively high densities and low temperatures, previously developed methods have been shown to be essential in improving equilibration by gradual insertion and deletion of these neutral sets of ions. We introduce an expanded ensemble approach which may be used in conjunction with these existing methods to further improve efficiency. Individual ions are inserted or deleted in one Monte Carlo trial rather than simultaneous insertion/deletion of neutral sets. We show how charge neutrality is maintained and show rigorous quantitative agreement between the conventional and the proposed expanded ensemble approaches, but with up to an order of magnitude increase in efficiency at high densities. The expanded ensemble approach is also more straightforward to implement than simultaneous insertion/deletion of neutral sets, and its implementation is demonstrated within open source software.

Construction of the interface potential from a series of canonical ensemble simulations.

We introduce a method to construct the interface potential from a series of molecular dynamics simulations conducted within the canonical ensemble. The interface potential provides the surface excess free energy associated with the growth of a fluid film from a surface. We collect the force that the fluid exerts on the surface (disjoining pressure) at a series of film thicknesses. These force data are then integrated to obtain the interface potential. "Spreading" and "drying" versions of the general approach are considered. The spreading approach focuses on the growth of a thin liquid film from a solid substrate in a mother vapor. The drying approach focuses on the growth of a thin vapor film on a solid substrate in a mother liquid. The methods provide a means to compute the contact angle of a fluid droplet in contact with the surface. The general method is applied to two model systems: (1) a monatomic Lennard-Jones fluid in contact with atomistically detailed face centered cubic (FCC) substrate and (2) TIP4P/2005 water in contact with a rigid silica surface. For the Lennard-Jones model system, we generate results with both the drying and spreading methods at various temperatures and substrate strengths. These results are compared to those from previous simulation studies. For the water system, the drying method is used to obtain wetting properties over a range of temperatures. The water system also highlights challenges associated with application of the spreading method within the framework pursued here.

Application of the interface potential approach for studying wetting behavior within a molecular dynamics framework.

We introduce a means to implement the interface potential approach for computing wetting properties within a molecular dynamics framework. The general approach provides a means to determine the contact angle of a liquid droplet on a solid substrate in a mother vapor. We present a framework for implementing "spreading" and "drying" versions of the method within an isothermal-isobaric ensemble. Two free energy methods are considered: cumulative integration of average force profile and multistate Bennett acceptance ratio. An umbrella sampling strategy is used to restrain volume fluctuations and to ensure adequate sampling of a broad volume range. We explore implementation of the approach with the GROningen MAchine for Chemical Simulations and the Large-scale Atomic/Molecular Massively Parallel Simulator. We test the accuracy and efficiency of the method with models consisting of a monoatomic Lennard-Jones fluid in the vicinity of a structureless or atomistically detailed substrate. Our results show that one can successfully generate the drying potential within the framework pursued here. The efficiency of the method is strongly dependent upon how one handles the dynamics of the two confining walls. These decisions impact the rate of volume fluctuations, and therefore, the quality of the volume distributions collected. Our efforts to implement the spreading method with molecular dynamics alone proved unsuccessful. The rate at which the configuration space of the vapor phase evolves is insufficient. We show how one can overcome this challenge by implementing a coupled molecular dynamics/Monte Carlo approach. Finally, we show how one can determine the variation in interfacial properties with temperature and substrate strength.

Effect of Carboxylic Acid on the Wetting Properties of a Model Water-Octane-Silica System.

Monte Carlo simulations are employed to determine the effects of acetic acid on the wetting properties of a model water-octane-silica system. We first compute the bulk liquid-vapor saturation properties of pure acetic acid and subsequently explore the bulk liquid-liquid saturation properties of the ternary water-octane-acid system. We introduce an expanded ensemble approach to compute the coexistence properties of the ternary system. An interface potential approach is then used to capture the evolution of the wetting properties of the water-octane-silica system upon the addition of acetic acid. We track the change in the octane-water interfacial tension and the contact angle of water droplet on a silica substrate in a mother octane fluid over a range of acetic acid activities. The structure of the fluid, including the partitioning of acetic acid within the interfacial system, is also considered at several state points. We observe that acetic acid has a strong tendency to adsorb at the octane-water interface, resulting in a reduction in the octane-water interfacial tension. The response of the contact angle is more sensitive to the temperature and the hydrophilicity of the silica substrate.

Using isothermal-isobaric Monte Carlo simulation to study the wetting behavior of model systems.

We introduce a molecular simulation method to compute the interfacial properties of model systems within the isothermal-isobaric ensemble. We use a free-energy-based approach in which Monte Carlo simulations are employed to obtain an interface potential associated with the growth of a fluid film from a solid substrate. The general method is implemented within "spreading" and "drying" frameworks. The interface potentials that emerge from these calculations provide direct access to spreading and drying coefficients. These macroscopic properties are then used to compute the liquid-vapor surface tension and the contact angle of a liquid droplet in contact with the substrate. The isothermal-isobaric ensemble provides a means to change the thickness of the fluid film adjacent to the substrate by modifying the volume of the simulation box. Molecular insertions and removals are not necessary. We introduce a framework for performing local volume change moves wherein one attempts to modify the density of a narrow region of the simulation box. We show that such moves improve the sampling efficiency of inhomogeneous systems. The approach is applied to a model system consisting of a monatomic Lennard-Jones fluid in the vicinity of a structureless substrate. Results are provided for direct spreading and drying interface potential calculations at several temperatures and substrate strengths. Expanded ensemble techniques are used to evaluate interfacial properties over a wide range of temperatures and substrate strengths. The results obtained using the isothermal-isobaric approach are compared with those previously obtained via a grand canonical approach.

Calculation of the Saturation Properties of a Model Octane-Water System Using Monte Carlo Simulation.

We use Monte Carlo simulation to compute the saturation properties of a model octane-water system. The system phase separates into water-rich liquid and vapor phases, octane-rich liquid and vapor phases, and water-rich liquid and octane-rich liquid phases at various conditions. We outline a strategy for determining the saturation properties of the mixture over a wide range of temperatures, pressures, and compositions. The approach begins with direct calculations that enable one to locate a single saturation point. A variety of expanded ensemble schemes are then used to trace saturation curves along paths of interest. We show how the overall strategy provides a means to construct pressure-composition diagrams at a fixed temperature and temperature-composition diagrams at a fixed pressure. In addition, we demonstrate how the approach is used to trace the liquid-liquid-vapor triple line over a wide range of temperatures. Simulation data are compared with experimental data, when available. Overall, our results show that the approach provides an efficient means to calculate the saturation properties of a binary system.

Multivariable extrapolation of grand canonical free energy landscapes.

We derive an approach for extrapolating the free energy landscape of multicomponent systems in the grand canonical ensemble, obtained from flat-histogram Monte Carlo simulations, from one set of temperature and chemical potentials to another. This is accomplished by expanding the landscape in a Taylor series at each value of the order parameter which defines its macrostate phase space. The coefficients in each Taylor polynomial are known exactly from fluctuation formulas, which may be computed by measuring the appropriate moments of extensive variables that fluctuate in this ensemble. Here we derive the expressions necessary to define these coefficients up to arbitrary order. In principle, this enables a single flat-histogram simulation to provide complete thermodynamic information over a broad range of temperatures and chemical potentials. Using this, we also show how to combine a small number of simulations, each performed at different conditions, in a thermodynamically consistent fashion to accurately compute properties at arbitrary temperatures and chemical potentials. This method may significantly increase the computational efficiency of biased grand canonical Monte Carlo simulations, especially for multicomponent mixtures. Although approximate, this approach is amenable to high-throughput and data-intensive investigations where it is preferable to have a large quantity of reasonably accurate simulation data, rather than a smaller amount with a higher accuracy.

Position-Dependent Dynamics Explain Pore-Averaged Diffusion in Strongly Attractive Adsorptive Systems.

Using molecular simulations, we investigate the relationship between the pore-averaged and position-dependent self-diffusivity of a fluid adsorbed in a strongly attractive pore as a function of loading. Previous work (Krekelberg, W. P.; Siderius, D. W.; Shen, V. K.; Truskett, T. M.; Errington, J. R. Connection between thermodynamics and dynamics of simple fluids in highly attractive pores. Langmuir 2013, 29, 14527-14535, doi: 10.1021/la4037327) established that pore-averaged self-diffusivity in the multilayer adsorption regime, where the fluid exhibits a dense film at the pore surface and a lower density interior pore region, is nearly constant as a function of loading. Here we show that this puzzling behavior can be understood in terms of how loading affects the fraction of particles that reside in the film and interior pore regions as well as their distinct dynamics. Specifically, the insensitivity of pore-averaged diffusivity to loading arises from the approximate cancellation of two factors: an increase in the fraction of particles in the higher diffusivity interior pore region with loading and a corresponding decrease in the particle diffusivity in that region. We also find that the position-dependent self-diffusivities scale with the position-dependent density. We present a model for predicting the pore-average self-diffusivity based on the position-dependent self-diffusivity, which captures the unusual characteristics of pore-averaged self-diffusivity in strongly attractive pores over several orders of magnitude.

Connection Between Thermodynamics and Dynamics of Simple Fluids in Pores: Impact of Fluid-Fluid Interaction Range and Fluid-Solid Interaction Strength.

Using molecular simulations, we investigate how the range of fluid-fluid (adsorbate-adsorbate) interactions and the strength of fluid-solid (adsorbate-adsorbent) interactions impact the strong connection between distinct adsorptive regimes and distinct self-diffusivity regimes reported in [Krekelberg, W. P.; Siderius, D. W.; Shen, V. K.; Truskett, T. M.; Errington, J. R. Langmuir2013, 29, 14527-14535]. Although increasing the fluid-fluid interaction range changes both the thermodynamics and the dynamic properties of adsorbed fluids, the previously reported connection between adsorptive filling regimes and self-diffusivity regimes remains. Increasing the fluid-fluid interaction range leads to enhanced layering and decreased self-diffusivity in the multilayer-formation regime but has little effect on the properties within film-formation and pore-filling regimes. We also find that weakly attractive adsorbents, which do not display distinct multilayer formation, are hard-sphere-like at super- and subcritical temperatures. In this case, the self-diffusivity of the confined and bulk fluid has a nearly identical scaling-relationship with effective density.

Temperature extrapolation of multicomponent grand canonical free energy landscapes.

We derive a method for extrapolating the grand canonical free energy landscape of a multicomponent fluid system from one temperature to another. Previously, we introduced this statistical mechanical framework for the case where kinetic energy contributions to the classical partition function were neglected for simplicity [N. A. Mahynski et al., J. Chem. Phys. 146, 074101 (2017)]. Here, we generalize the derivation to admit these contributions in order to explicitly illustrate the differences that result. Specifically, we show how factoring out kinetic energy effects a priori, in order to consider only the configurational partition function, leads to simpler mathematical expressions that tend to produce more accurate extrapolations than when these effects are included. We demonstrate this by comparing and contrasting these two approaches for the simple cases of an ideal gas and a non-ideal, square-well fluid.

Predicting low-temperature free energy landscapes with flat-histogram Monte Carlo methods.

We present a method for predicting the free energy landscape of fluids at low temperatures from flat-histogram grand canonical Monte Carlo simulations performed at higher ones. We illustrate our approach for both pure and multicomponent systems using two different sampling methods as a demonstration. This allows us to predict the thermodynamic behavior of systems which undergo both first order and continuous phase transitions upon cooling using simulations performed only at higher temperatures. After surveying a variety of different systems, we identify a range of temperature differences over which the extrapolation of high temperature simulations tends to quantitatively predict the thermodynamic properties of fluids at lower ones. Beyond this range, extrapolation still provides a reasonably well-informed estimate of the free energy landscape; this prediction then requires less computational effort to refine with an additional simulation at the desired temperature than reconstruction of the surface without any initial estimate. In either case, this method significantly increases the computational efficiency of these flat-histogram methods when investigating thermodynamic properties of fluids over a wide range of temperatures. For example, we demonstrate how a binary fluid phase diagram may be quantitatively predicted for many temperatures using only information obtained from a single supercritical state.

Understanding the influence of Coulomb and dispersion interactions on the wetting behavior of ionic liquids.

We study the role of dispersion and electrostatic interactions in the wetting behavior of ionic liquids on non-ionic solid substrates. We consider a simple model of an ionic liquid consisting of spherical ions that interact via Lennard-Jones and Coulomb potentials. Bulk and interfacial properties are computed for five fluids distinguished by the strength of the electrostatic interaction relative to the dispersion interaction. We employ Monte Carlo simulations and an interface-potential-based approach to calculate the liquid-vapor and substrate-fluid interfacial properties. Surface tensions for each fluid are evaluated over a range of temperatures that spans from a reduced temperature of approximately 0.6 to the critical point. Contact angles are calculated at select temperatures over a range of substrate-fluid interaction strengths that spans from the near-drying regime to the wetting regime. We observe that an increase in the relative strength of Coulombic interactions between ions leads to increasing deviation from Guggenheim's corresponding states theory. We show how this deviation is related to lower values of liquid-vapor excess entropies observed for strongly ionic fluids. Our results show that the qualitative nature of wetting behavior is significantly influenced by the competition between dispersion and electrostatic interactions. We discuss the influence of electrostatic interactions on the nature of wetting and drying transitions and corresponding states like behavior observed for contact angles. For all of the fluids studied, we observe a relatively narrow range of substrate-fluid interaction strengths wherein the contact angle is nearly independent of temperature. The influence of the ionic nature of the fluid on the temperature dependence of contact angle is also discussed.

Saturation properties of 1-alkyl-3-methylimidazolium based ionic liquids.

We study the liquid-vapor saturation properties of room temperature ionic liquids (RTILs) belonging to the homologous series 1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([Cnmim][NTf2]) using Monte Carlo simulation. We examine the effect of temperature and cation alkyl chain length n on the saturated densities, vapor pressures, and enthalpies of vaporization. These properties are explicitly calculated for temperatures spanning from 280 to 1000 K for RTILs with n = 2, 4, 6, 8, 10, and 12. We also explore how the identity of the anion influences saturation properties. Specifically, we compare results for [C(4)mim][NTf2] with those for 1-butyl-3-methylimidazolium tetrafluoroborate ([C(4)mim][BF4]) and 1-butyl-3-methylimidazolium hexafluorophosphate ([C(4)mim][PF6]). Simulations are completed with a recently developed realistic united-atom force field. A combination of direct grand canonical and isothermal-isobaric temperature expanded ensemble simulations are used to construct phase diagrams. Our results are compared with experimental data and Gibbs ensemble simulation data. Overall, we find good agreement between our results and those measured experimentally. We find that the vapor pressures and enthalpies of vaporization show a strong dependence on the size of the alkyl chain at low temperatures, whereas no particular trend is observed at high temperatures. Finally, we also discuss the effect of temperature on liquid phase nanodomains observed in RTILs with large hydrophobic groups. We do not observe a drastic change in liquid phase structure upon variation of the temperature, which suggests there is not a sharp phase transition between a nanostructured and homogeneous liquid, as has been suggested in earlier studies.

Connection between thermodynamics and dynamics of simple fluids in highly attractive pores.

Using molecular simulations, we investigate the structural and diffusive dynamics properties of a model fluid in highly absorptive cylindrical pores. At subcritical temperatures, self-diffusivity displays three distinct regimes as a function of average pore density ρ: (1) a decrease in self-diffusivity with increasing ρ at low ρ, (2) constant self-diffusivity with respect to varying ρ at moderate density, and (3) a decrease in self-diffusivity with increasing ρ at high density. These regimes are closely linked to the thermodynamic properties of the fluid in the pore, specifically, the adsorption isotherm, isosteric heat of adsorption, and the density profile. We show that these three diffusivity regimes qualitatively correspond to three distinct adsorption regimes: monolayer formation, multilayer adsorption, and pore filling, respectively. In addition, we find that the self-diffusivity is a universal function of the local film density in the monolayer formation regime at subcritical temperatures. The results of this work suggest a potential means to estimate the self-diffusivity over a broad pressure range using a limited number of experiments.

Communication: Phase behavior of materials with isotropic interactions designed by inverse strategies to favor diamond and simple cubic lattice ground states.

We use molecular simulation to construct equilibrium phase diagrams for two recently introduced model materials with isotropic, soft-repulsive pair interactions designed to favor diamond and simple cubic lattice ground states, respectively, over a wide range of densities [Jain et al., Soft Matter 9, 3866 (2013)]. We employ free energy based Monte Carlo simulation techniques to precisely trace the inter-crystal and fluid-crystal coexistence curves. We find that both model materials display rich polymorphic phase behavior featuring stable crystals corresponding to the target ground-state structures, as well as a variety of other crystalline (e.g., hexagonal and body-centered cubic) phases and multiple reentrant melting transitions.

Impact of small-scale geometric roughness on wetting behavior.

We examine the extent to which small-scale geometric substrate roughness influences the wetting behavior of fluids at solid surfaces. Molecular simulation is used to construct roughness wetting diagrams wherein the progression of the contact angle is traced from the Cassie to Wenzel to impregnation regime with increasing substrate strength for a collection of systems with rectangularly shaped grooves. We focus on the evolution of these diagrams as the length scale of the substrate features approaches the size of a fluid molecule. When considering a series of wetting diagrams for substrates with fixed shape and variable feature periodicity, we find that the diagrams progressively shift away from a common curve as the substrate features become smaller than approximately 10 fluid diameters. It is at this length scale that the macroscopic models of Cassie and Wenzel become unreliable. Deviations from the macroscopic models are attributed to the manner in which the effective substrate-fluid interaction strength evolves with periodicity and the important role that confinement effects play for substrates with small periodicities.

Understanding wetting of immiscible liquids near a solid surface using molecular simulation.

We introduce Monte Carlo simulation methods for determining interfacial properties of fluids that exhibit bulk liquid-liquid immiscibility. 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, is utilized to deduce the wetting characteristics of these systems. We present a framework for implementing this general method within both the grand canonical and semigrand isobaric-isothermal ensembles. Tracking the evolution of interfacial properties along various thermodynamic paths is also examined. This task is accomplished by implementing variants of the expanded ensemble technique, which enables one to obtain components of the interface potential along a path of interest. We also discuss how these concepts are employed to calculate bulk liquid-liquid coexistence properties in an efficient manner. The computational strategies introduced here are applied to three model Lennard-Jones systems. For each system, we compile the evolution of the liquid-liquid surface tension and contact angle with temperature or pressure. For one of the model systems we compare our results with literature data. We also examine how interfacial properties evolve upon variation of the relative affinity of the fluid components for the substrate. Overall, we find that the approach pursued here is generally applicable and provides an efficient and precise means to calculate the bulk and interfacial properties of fluids that exhibit liquid-liquid immiscibility.

Using Monte Carlo simulation to compute liquid-vapor saturation properties of ionic liquids.

We discuss Monte Carlo (MC) simulation methods for calculating liquid-vapor saturation properties of ionic liquids. We first describe how various simulation tools, including reservoir grand canonical MC, growth-expanded ensemble MC, distance-biasing, and aggregation-volume-biasing, are used to address challenges commonly encountered in simulating realistic models of ionic liquids. We then indicate how these techniques are combined with histogram-based schemes for determining saturation properties. Both direct methods, which enable one to locate saturation points at a given temperature, and temperature expanded ensemble methods, which provide a means to trace saturation lines to low temperature, are discussed. We study the liquid-vapor phase behavior of the restricted primitive model (RPM) and a realistic model for 1,3-dimethylimidazolium tetrafluoroborate ([C1mim][BF4]). Results are presented to show the dependence of saturation properties of the RPM and [C1mim][BF4] on the size of the simulation box and the boundary condition used for the Ewald summation. For [C1mim][BF4] we also demonstrate the ability of our strategy to sample ion clusters that form in the vapor phase. Finally, we provide the liquid-vapor saturation properties of these models over a wide range of temperature. Overall, we observe that the choice of system size and boundary condition have a non-negligible effect on the calculated properties, especially at high temperature. Also, we find that the combination of grand canonical MC simulation and isothermal-isobaric temperature expanded ensemble MC simulation provides a computationally efficient means to calculate liquid-vapor saturation properties of ionic liquids.