Coexistence diagrams of mixtures by molecular simulation

Modelling Sorption and Transport of Gases in Polymeric Membranes across Different Scales: A Review

Professor Giulio C. Sarti has provided outstanding contributions to the modelling of fluid sorption and transport in polymeric materials, with a special eye on industrial applications such as membrane separation, due to his Chemical Engineering background. He was the co-creator of innovative theories such as the Non-Equilibrium Theory for Glassy Polymers (NET-GP), a flexible tool to estimate the solubility of pure and mixed fluids in a wide range of polymers, and of the Standard Transport Model (STM) for estimating membrane permeability and selectivity. In this review, inspired by his rigorous and original approach to representing membrane fundamentals, we provide an overview of the most significant and up-to-date modeling tools available to estimate the main properties governing polymeric membranes in fluid separation, namely solubility and diffusivity. The paper is not meant to be comprehensive, but it focuses on those contributions that are most relevant or that show the potential to be relevant in the future. We do not restrict our view to the field of macroscopic modelling, which was the main playground of professor Sarti, but also devote our attention to Molecular and Multiscale Hierarchical Modeling. This work proposes a critical evaluation of the different approaches considered, along with their limitations and potentiality.

On the calculation of free energies over Hamiltonian and order parameters via perturbation and thermodynamic integration

In this work, complementary formulas are presented to compute free-energy differences via perturbation (FEP) methods and thermodynamic integration (TI). These formulas are derived by selecting only the most statistically significant data from the information extractable from the simulated points involved. On the one hand, commonly used FEP techniques based on overlap sampling leverage the full information contained in the overlapping macrostate probability distributions. On the other hand, conventional TI methods only use information on the first moments of those distributions, as embodied by the first derivatives of the free energy. Since the accuracy of simulation data degrades considerably for high-order moments (for FEP) or free-energy derivatives (for TI), it is proposed to consider, consistently for both methods, data up to second-order moments/derivatives. This provides a compromise between the limiting strategies embodied by common FEP and TI and leads to simple, optimized expressions to evaluate free-energy differences. The proposed formulas are validated with an analytically solvable harmonic Hamiltonian (for assessing systematic errors), an atomistic system (for computing the potential of mean force with coordinate-dependent order parameters), and a binary-component coarse-grained model (for tracing a solid-liquid phase diagram in an ensemble sampled through alchemical transformations). It is shown that the proposed FEP and TI formulas are straightforward to implement, perform similarly well, and allow robust estimation of free-energy differences even when the spacing of successive points does not guarantee them to have proper overlapping in phase space.

Direct evaluation of the phase diagrams of dense multicomponent plasmas by integration of the Clapeyron equations

Accurate phase diagrams of multicomponent plasmas are required for the modeling of dense stellar plasmas, such as those found in the cores of white dwarf stars and the crusts of neutron stars. Those phase diagrams have been computed using a variety of standard techniques, which suffer from physical and computational limitations. Here we present an efficient and accurate method that overcomes the drawbacks of previously used approaches. In particular, finite-size effects are avoided as each phase is calculated separately; the plasma electrons and volume changes are explicitly taken into account; and arbitrary analytic fits to simulation data as well as particle insertions are avoided. Furthermore, no simulations at "uninteresting" state conditions, i.e., away from the phase coexistence curves, are required, which improves the efficiency of the technique. The method consists of an adaptation of the so-called Gibbs-Duhem integration approach to electron-ion plasmas, where the coexistence curve is determined by direct numerical integration of its underlying Clapeyron equation. The thermodynamics properties of the coexisting phases are evaluated separately using Monte Carlo simulations in the isobaric semigrand canonical ensemble (NPTΔμ). We describe this Monte Carlo-based Clapeyron integration method, including its basic physical and numerical principles, our extension to electron-ion plasmas, and our numerical implementation. We illustrate its applicability and benefits with the calculation of the melting curve of dense carbon-oxygen plasmas under conditions relevant for the cores of white dwarf stars and provide analytic fits to implement this new melting curve in white dwarf models. While this work focuses on the liquid-solid phase boundary of dense two-component plasmas, a wider range of physical systems and phase boundaries are within the scope of the Clapeyron integration method, which had until now only been applied to simple model systems of neutral particles.

Efficient determination of solid-state phase equilibrium with the multicell Monte Carlo method

Building on our previously introduced multicell Monte Carlo (MC)^{2} method for modeling phase coexistence, this paper provides important improvements for efficient determination of phase equilibria in solids. The (MC)^{2} method uses multiple cells, representing possible phases. Mass transfer between cells is modeled virtually by solving the mass balance equation after the composition of each cell is changed arbitrarily. However, searching for the minimum free energy during this process poses a practical problem. The solution to the mass balance equation is not unique away from equilibrium, and consequently the algorithm is in risk of getting trapped in nonequilibrium solutions. Therefore, a proper stopping condition for (MC)^{2} is currently lacking. In this work, we introduce a consistency check via a predictor-corrector algorithm to penalize solutions that do not satisfy a necessary condition for equivalence of chemical potentials and steer the system toward finding equilibrium. The most general acceptance criteria for (MC)^{2} is derived starting from the isothermal-isobaric Gibbs ensemble for mixtures. Using this ensemble, translational MC moves are added to include vibrational excitations as well as volume MC moves to ensure the condition of constant pressure and temperature entirely with a MC approach, without relying on any other method for relaxation of these degrees of freedom. As a proof of concept the method is applied to two binary alloys with miscibility gaps and a model quaternary alloy, using classical interatomic potentials.

On flexible force fields for metal-organic frameworks: Recent developments and future prospects

Classical force field simulations can be used to study structural, diffusion, and adsorption properties of metal-organic frameworks (MOFs). To account for the dynamic behavior of the material, parameterization schemes have been developed to derive force constants and the associated reference values by fitting on ab initio energies, vibrational frequencies, and elastic constants. Here, we review recent developments in flexible force field models for MOFs. Existing flexible force field models are generally able to reproduce the majority of experimentally observed structural and dynamic properties of MOFs. The lack of efficient sampling schemes for capturing stimuli-driven phase transitions, however, currently limits the full predictive potential of existing flexible force fields from being realized. This article is categorized under: Structure and Mechanism > Computational Materials ScienceMolecular and Statistical Mechanics > Molecular Mechanics.

Influence of patch-size variability on the crystallization of tetrahedral patchy particles

The understanding of disorder effects on crystallization is of fundamental and technological importance. It is well established by both theory and experiment that particle-size polydispersity hinders crystallization for isotropically interacting particles. Here, we address the effects of patch variability in a model for tetrahedral colloids, where polydispersity is introduced independently on the size, position, and strength of the attractive patches. Our simulations indicate that, unlike particle-size polydispersity, angular polydispersity has a minor impact on the crystallization properties of tetrahedral colloidal particles. Particles with angular polydispersity well within current experimental possibilities fully retain their crystallization properties, a result which should encourage the realization of colloidal crystals in experiment.

Coexistence of fluid and crystalline phases of proteins in photosynthetic membranes

Photosystem II (PSII) and its associated light-harvesting complex II (LHCII) are highly concentrated in the stacked grana regions of photosynthetic thylakoid membranes. PSII-LHCII supercomplexes can be arranged in disordered packings, ordered arrays, or mixtures thereof. The physical driving forces underlying array formation are unknown, complicating attempts to determine a possible functional role for arrays in regulating light harvesting or energy conversion efficiency. Here, we introduce a coarse-grained model of protein interactions in coupled photosynthetic membranes, focusing on just two particle types that feature simple shapes and potential energies motivated by structural studies. Reporting on computer simulations of the model's equilibrium fluctuations, we demonstrate its success in reproducing diverse structural features observed in experiments, including extended PSII-LHCII arrays. Free energy calculations reveal that the appearance of arrays marks a phase transition from the disordered fluid state to a system-spanning crystal. The predicted region of fluid-crystal coexistence is broad, encompassing much of the physiologically relevant parameter regime; we propose experiments that could test this prediction. Our results suggest that grana membranes lie at or near phase coexistence, conferring significant structural and functional flexibility to this densely packed membrane protein system.

Thermodynamic methods and models to study flexible metal-organic frameworks

Much attention has recently been focused on a fascinating subclass of metal-organic frameworks that behave in a remarkable stimuli-responsive fashion. These soft porous crystals feature dynamic crystalline frameworks displaying reversible, large-amplitude structural deformations under external physical constraints such as temperature, electric field or gas exposure. The number of reported syntheses of such materials is rapidly growing and they are promising for practical applications, such as gas capture, purification and fluid separation. Herein, we summarize the recently developed thermodynamic tools that can help understand the process of fluid adsorption and fluid mixture coadsorption in these flexible nanoporous materials. These tools, which include both molecular simulation methods and analytical models, can help rationalize experimental results and predict adsorption properties over a wide range of thermodynamic conditions. A particular focus is given on how these methods can guide the experimental exploration of a large number of materials and working conditions (temperature, pressure, composition) to help design efficient processes relying on fluid adsorption in soft porous crystals.

Virial coefficients of Lennard-Jones mixtures

We report results of calculations of the second through sixth virial coefficients for four prototype Lennard-Jones (LJ) mixtures that have been the subject of previous studies in the literature. Values are reported for temperatures ranging from T=0.6 to T=10.0, where here the temperature is given units of the LJ energy parameter of one of the components. Thermodynamic stability of the mixtures is studied using the virial equation of state (VEOS) with the calculated coefficients, with particular focus on characterizing the vapor-liquid critical behavior of the mixtures. For three of the mixtures, vapor-liquid coexistence and critical data are available for comparison at only one temperature, while for the fourth we can compare to a critical line. We find that the VEOS provides a useful indication of the presence and location of critical behavior, although in some situations we find need to consider "near-miss" critical behavior, where the classical conditions of criticality are nearly but not exactly satisfied.

Improvement in molecule exchange efficiency in Gibbs ensemble Monte Carlo: development and implementation of the continuous fractional component move

The continuous fractional component Monte Carlo (CFC MC) move (J Chem Theory Comput, 2007, 3, 1451) is extended to the Gibbs ensemble. The algorithm is validated against conventional simulations for the Lennard Jones fluid and a flexible water model. The method is also used to compute the vapor-liquid coexistence densities of a model for SO(2). The CFC molecule exchange move relies on the gradual insertion and deletion of molecules in conjunction with a self-adapting bias potential. As a result, the method does not require the formation of spontaneous voids in the dense fluid phase to be successful, leading to molecule exchange acceptance probabilities that are nearly independent of temperature. For example, over 1% of the vapor-liquid molecule exchange moves are successful for water at 280 K, whereas advanced rotational and configurational bias methods have success rates of less than 0.03%. The CFC move can be combined with other Monte Carlo moves to enable efficient simulation of dense strongly associating fluids that are to this point infeasible to model with standard methods.

Optimized expanded ensembles for simulations involving molecular insertions and deletions. II. Open systems

In the Grand Canonical, osmotic, and Gibbs ensembles, chemical potential equilibrium is attained via transfers of molecules between the system and either a reservoir or another subsystem. In this work, the expanded ensemble (EXE) methods described in part I [F. A. Escobedo and F. J. Martinez-Veracoechea, J. Chem. Phys. 127, 174103 (2007)] of this series are extended to these ensembles to overcome the difficulties associated with implementing such whole-molecule transfers. In EXE, such moves occur via a target molecule that undergoes transitions through a number of intermediate coupling states. To minimize the tunneling time between the fully coupled and fully decoupled states, the intermediate states could be either: (i) sampled with an optimal frequency distribution (the sampling problem) or (ii) selected with an optimal spacing distribution (staging problem). The sampling issue is addressed by determining the biasing weights that would allow generating an optimal ensemble; discretized versions of this algorithm (well suited for small number of coupling stages) are also presented. The staging problem is addressed by selecting the intermediate stages in such a way that a flat histogram is the optimized ensemble. The validity of the advocated methods is demonstrated by their application to two model problems, the solvation of large hard spheres into a fluid of small and large spheres, and the vapor-liquid equilibrium of a chain system.

Isomolar-semigrand ensemble molecular dynamics: Application to vapor-liquid equilibrium of the mixture methane/ethane

The isomolar-semigrand ensemble molecular dynamics (iSGMD) method is applied to the simulation of the binary system methane/ethane. The vapor-liquid equilibrium properties of this system at a temperature of 192.37 K are computed using the Gibbs-Duhem integration method. The iSGMD method, which resembles conventional hybrid Monte Carlo (MC) but is applicable to phase equilibrium calculations, is designed to overcome the difficulties associated with performing standard Monte Carlo-type particle transformations in liquid systems that are very dense and/or are comprised of complex molecules with many intramolecular degrees of freedom. This work shows that particle transformations using the iSGMD method for the simple system methane/ethane are at least 25 times more successful than standard MC-type transformations. The P-x-y curve for the system methane/ethane at 192.37 K computed using iSGMD simulations agrees very well with the experimental P-x-y curve as well as results of a previous MC study.

Three-body interactions and solid-liquid phase equilibria: application of a molecular dynamics algorithm

The effect of three-body interactions on the solid-liquid phase boundaries of argon, krypton, and xenon is investigated via a novel technique that combines both nonequilibrium and equilibrium molecular dynamics. The simulations involve the evaluation of two- and three-body forces using accurate two-body and three-body intermolecular potentials. The effect of three-body interactions is to substantially increase the coexistence pressure and to lower the densities of liquid and solid phases. Comparison with experiment indicates that three-body interactions are required to accurately determine the total pressure. In contrast to vapor-liquid phase equilibria, the relative contribution of three-body interactions to the freezing pressure exceeds the contribution of two-body interactions at all temperatures.

Effect of three-body interactions on the vapor-liquid phase equilibria of binary fluid mixtures

Gibbs-Duhem Monte Carlo simulations are reported for the vapor-liquid phase coexistence of binary argon+krypton mixtures at different temperatures. The calculations employ accurate two-body potentials in addition to contributions from three-body dispersion interactions resulting from third-order triple-dipole interactions. A comparison is made with experiment that illustrates the role of three-body interactions on the phase envelope. In all cases the simulations represent genuine predictions with input parameters obtained independently from sources other than phase equilibria data. Two-body interactions alone are insufficient to adequately describe vapor-liquid coexistence. In contrast, the addition of three-body interactions results in very good agreement with experiment. In addition to the exact calculation of three-body interactions, calculations are reported with an approximate formula for three-body interactions, which also yields good results.

Direct determination of liquid phase coexistence by Monte Carlo simulations

A formalism to determine coexistence points by means of Monte Carlo simulations is presented. The general idea of the method is to perform a simulation simultaneously in several unconnected boxes which can exchange particles. At equilibrium, most of the boxes will be occupied by a homogeneous phase. The compositions of these boxes yield coexisting points on the binodal. However, since the overall composition is fixed, at least one of the boxes will contain an interface. We show that this does not affect the results, provided that the interface has no net curvature. We coin the name "Helmholtz-ensemble method," because the method is related to the well-known Gibbs-ensemble method, but the volume of the boxes is constant. Since the box volumes are constant, we expect that this method will be particularly useful for lattice models. The accuracy of the Helmholtz-ensemble method is benchmarked against known coexistence curves of the three-dimensional Ising model with excellent results.

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

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

A general perturbation approach for equation of state development: Applications to simple fluids, ab initio potentials, and fullerenes

A new perturbation scheme based on the Barker-Henderson perturbation theory [J. Chem. Phys. 47, 4714 (1967)] is proposed to predict the thermodynamic properties of spherical molecules. Accurate predictions of second virial coefficients and vapor-liquid coexistence properties are obtained for a large variety of potential functions (square well, Yukawa, Sutherland, Lennard-Jones, Buckingham, Girifalco). New Gibbs ensemble Monte Carlo simulations of the generalized exp-m Buckingham potential are reported. An extension of the perturbation approach to mixtures is proposed, and excellent predictions of vapor-liquid equilibria are obtained for Lennard-Jones mixtures. The perturbation scheme can be applied to complex potential functions fitted to ab initio data to predict the properties of real molecules such as neon. The new approach can also be used as an auxiliary tool in molecular simulation studies, to efficiently optimize an intermolecular potential on macroscopic properties or match force fields based on different potential functions.

Computing the starting state for Gibbs-Duhem integration

Gibbs-Duhem integration implies the numerical integration of a Clapeyron equation. To start the numerical integration, an initial coexistence point and a corresponding initial slope of the Clapeyron equation are needed. In order to apply Gibbs-Duhem integration to all kinds of systems at diverse physical conditions, one has to investigate and assess the available methods that can be used to compute these initial values. This publication focuses on vapor-liquid equilibria in binary mixtures comprising chain molecules. The initial coexistence point is either computed with the NVbeta Gibbs ensemble or with the Npbeta+test molecule method with overlapping distributions, which is introduced in this publication. Although computationally demanding, the Npbeta+test molecule method with overlapping distributions is applicable at conditions where the NVbeta Gibbs ensemble fails. We investigated three methods that can be employed to compute the initial slope of the Clapeyron equation. The Widom method and the overlapping-distributions difference method provide correct values for the initial slope. The difference method does only provide the correct answer in special cases. The possibility to judge the reliability of the results makes the overlapping-distributions difference method the safest route to the initial slope. Gibbs-Duhem integration requires the frequent computation of the slope of the Clapeyron equation. This slope depends on ensemble averages of the composition. A new bias method for efficient sampling of the composition in a semigrand-canonical simulation of chain molecules is presented. This bias method considerably enhances the composition sampling in systems comprising chain molecules of different sizes.

An advanced Gibbs-Duhem integration method: Theory and applications

The conventional Gibbs-Duhem integration method is very convenient for the prediction of phase equilibria of both pure components and mixtures. However, it turns out to be inefficient. The method requires a number of lengthy simulations to predict the state conditions at which phase coexistence occurs. This number is not known from the outset of the numerical integration process. Furthermore, the molecular configurations generated during the simulations are merely used to predict the coexistence condition and not the liquid- and vapor-phase densities and mole fractions at coexistence. In this publication, an advanced Gibbs-Duhem integration method is presented that overcomes above-mentioned disadvantage and inefficiency. The advanced method is a combination of Gibbs-Duhem integration and multiple-histogram reweighting. Application of multiple-histogram reweighting enables the substitution of the unknown number of simulations by a fixed and predetermined number. The advanced method has a retroactive nature; a current simulation improves the predictions of previously computed coexistence points as well. The advanced Gibbs-Duhem integration method has been applied for the prediction of vapor-liquid equilibria of a number of binary mixtures. The method turned out to be very convenient, much faster than the conventional method, and provided smooth simulation results. As the employed force fields perfectly predict pure-component vapor-liquid equilibria, the binary simulations were very well suitable for testing the performance of different sets of combining rules. Employing Lorentz-Hudson-McCoubrey combining rules for interactions between unlike molecules, as opposed to Lorentz-Berthelot combining rules for all interactions, considerably improved the agreement between experimental and simulated data.

Determination of surface tension in binary mixtures using transition-matrix Monte Carlo

We present a methodology based on grand-canonical transition-matrix Monte Carlo and finite-size scaling analysis to calculate surface tensions in binary mixtures. In particular, mixture transition-matrix Monte Carlo is first used to calculate apparent, system-size-dependent free-energy barriers separating coexisting fluid phases. Finite-size scaling is then used to extrapolate these values to the infinitely large system limit to determine the true thermodynamic surface tension. A key distinction of the methodology is that it yields the entire isothermal surface-tension curve for a binary mixture in a relatively small number of simulations. We demonstrate the utility of the method by calculating surface-tension curves for three binary Lennard-Jones mixtures. While we have only examined the surface tension of simple fluids in this work, the method is general and can be extended to molecular fluids as well as to determine interfacial tensions of liquid-liquid interfaces.

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

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