Gibbs Ensemble Monte Carlo Simulations Research Articles

Optimising the parameters of the Gibbs Ensemble Monte Carlo simulation of phase separation: the role of multiple relaxation times

ABSTRACTWe address the problem of how to choose the parameters of a Gibbs Ensemble Monte Carlo (GEMC) simulation of phase separation optimally through a systematic study of the autocorrelation functions of density fluctuation in a particular example. In this model the particles interact via a square well potential. Values of the temperature as well as the range of the square well potential are varied. We find that the normalised autocorrelation function is described accurately over a very large time scale as a linear combination of an exponential function with a time scale τ2 and a stretched exponential function with a time scale τ1 and an exponent α. Dependence of α, τ1 and τ2 on the temperature, the parameters of the GEMC algorithm and the range of the square well potential is investigated and interpreted. We use the insight thus obtained and the generality of some aspects of these results to construct the route to optimisation of the parameters of a GEMC simulation in a somewhat more generic situation involving a single-component phase separation.

Molecular Mechanism for Azeotrope Formation in Ethanol/Benzene Binary Mixtures through Gibbs Ensemble Monte Carlo Simulation.

Azeotropes have been studied for decades due to the challenges they impose on separation processes but fundamental understanding at the molecular level remains limited. Although molecular simulation has demonstrated its capability of predicting mixture vapor-liquid equilibrium (VLE) behaviors, including azeotropes, its potential for mechanistic investigation has not been fully exploited. In this study, we use the united atom transferable potentials for phase equilibria (TraPPE-UA) force field to model the ethanol/benzene mixture, which displays a positive azeotrope. Gibbs ensemble Monte Carlo (GEMC) simulation is performed to predict the VLE phase diagram, including an azeotrope point. The results accurately agree with experimental measurements. We argue that the molecular mechanism of azeotrope formation cannot be fully understood by studying the mixture liquid-state stability at the azeotrope point alone. Rather, azeotrope occurrence is only a reflection of the changing relative volatility between the two components over a much wider composition range. A thermodynamic criterion is thus proposed on the basis of the comparison of partial excess Gibbs energy between the components. In the ethanol/benzene system, molecular energetics shows that with increasing ethanol mole fraction, its volatility initially decreases but later plateaus, while benzene volatility is initially nearly constant and only starts to decrease when its mole fraction is low. Analysis of the mixture liquid structure, including a detailed investigation of ethanol hydrogen-bonding configurations at different composition levels, reveals the underlying molecular mechanism for the changing volatilities responsible for the azeotrope.

Two-Phase Equilibrium Conditions in Nanopores.

It is known that thermodynamic properties of a system change upon confinement. To know how, is important for modelling of porous media. We propose to use Hill’s systematic thermodynamic analysis of confined systems to describe two-phase equilibrium in a nanopore. The integral pressure, as defined by the compression energy of a small volume, is then central. We show that the integral pressure is constant along a slit pore with a liquid and vapor in equilibrium, when Young and Young–Laplace’s laws apply. The integral pressure of a bulk fluid in a slit pore at mechanical equilibrium can be understood as the average tangential pressure inside the pore. The pressure at mechanical equilibrium, now named differential pressure, is the average of the trace of the mechanical pressure tensor divided by three as before. Using molecular dynamics simulations, we computed the integral and differential pressures, and p, respectively, analysing the data with a growing-core methodology. The value of the bulk pressure was confirmed by Gibbs ensemble Monte Carlo simulations. The pressure difference times the volume, V, is the subdivision potential of Hill, . The combined simulation results confirm that the integral pressure is constant along the pore, and that scales with the inverse pore width. This scaling law will be useful for prediction of thermodynamic properties of confined systems in more complicated geometries.

Molecular modeling of vapor-liquid equilibrium properties of HFC-161 and its mixture HFC-161+HFO-1234yf

Abstract Molecular simulation is an important method to acquire reliable thermophysical properties, but there are few force fields for alternative refrigerants. In this paper, an all-atom force field for ethyl fluoride (HFC-161) has been developed in the framework of the Amber force field. This force field for HFC-161 was derived from a combination of ab initio calculations, Gibbs ensemble Monte Carlo simulations and comparisons to experimental data. And then the force field was applied to the calculation of thermophysical properties of pure HFC-161. It was found that over the temperature range of 263.15 K to 363.15 K the average absolute relative deviations of vapor pressure, saturated densities in liquid and vapor phases were 1.37%, 3.87% and 1.86%, respectively. The deviations in critical temperature, critical pressure and critical density were 0.40%, 1.86% and 1.47%, respectively. This force field was also used to calculate the vapor-liquid equilibrium properties of the mixture of HFC-161 + HFO-1234yf. In this work, the vapor-liquid equilibrium properties of pure HFC-161 and its binary mixture were calculated via the Gibbs ensemble Monte Carlo simulations. Calculated results showed a good agreement with available experimental data. Since the calculated results were in good agreement with the experimental results, the predictive ability of the proposed force field model can be verified.

Gibbs Ensemble Monte Carlo Simulation of Fluids in Confinement: Relation between the Differential and Integral Pressures.

The accurate description of the behavior of fluids in nanoporous materials is of great importance for numerous industrial applications. Recently, a new approach was reported to calculate the pressure of nanoconfined fluids. In this approach, two different pressures are defined to take into account the smallness of the system: the so-called differential and the integral pressures. Here, the effect of several factors contributing to the confinement of fluids in nanopores are investigated using the definitions of the differential and integral pressures. Monte Carlo (MC) simulations are performed in a variation of the Gibbs ensemble to study the effect of the pore geometry, fluid-wall interactions, and differential pressure of the bulk fluid phase. It is shown that the differential and integral pressure are different for small pores and become equal as the pore size increases. The ratio of the driving forces for mass transport in the bulk and in the confined fluid is also studied. It is found that, for small pore sizes (i.e., <), the ratio of the two driving forces considerably deviates from 1.

Diagram of vapor‐liquid equilibria for n‐pentane using hybrid Gibbs ensemble Monte Carlo simulation

AbstractIn this paper we report the simulation results using the 5‐site Lennard‐Jones potential function obtained from ab initio interaction energy for n‐pentane at level of theory CCSD(T)/aug‐cc‐pVTZ. The hybrid Gibbs ensemble Monte Carlo (HGEMC) simulation was performed under the same initial conditions for boxes 1 and 2 with 60 molecules and a density of 0.130 g /cm3. The HGEMC simulation was implemented with 106 epochs and value of Cutoff_Radius 10Å. The simulation temperature was varied from 100 to 450 K. The thermodynamic results of the vapor‐liquid equilibria derived from HGEMC‐NVT simulation are consistent with experimental data and other calculation methods.

The effect of surface roughness on the phase behavior of colloidal particles.

Shape anisotropy of colloidal particles can give rise to complex intermolecular interactions that determine particle packing and phase behavior. The vapor-liquid coexistence curves of attractive rough particles display a shift when compared to attractive smooth spherical particles. We use Integral Equation Theory (IET) to determine the vapor-liquid spinodal phase diagram of smooth and rough colloidal particles interacting through square-well attraction. Additionally, we use Gibbs Ensemble Monte Carlo (GEMC) simulations to locate their vapor-liquid coexistence curves. We model a rough colloidal particle as a spherical core with small beads embedded on its surface. The critical point of smooth spherical particle systems predicted by theory and simulations is in quantitative agreement. An increase in surface roughness due to an increase in either the number of beads or the diameter of the beads has a modest effect on the local structure of the system in the supercritical region. In contrast, increasing surface roughness consistently shifts the vapor-liquid coexistence curves to higher temperatures. The critical temperature is found to be a quadratic function of the number of beads. At a fixed bead size and number of beads, the critical temperature does not vary with the arrangement of beads on the core. Both IET and GEMC simulations predict that unlike critical temperatures, critical packing fractions vary non-monotonically with surface roughness. We find that the feasibility and accuracy of the integral equation theory depend sensitively on the chosen closure combination.

Thermodynamic properties of the 3D Lennard-Jones/spline model

The Lennard-Jones (LJ) spline potential is a truncated LJ potential so that both the pair potential and the force continuously approach zero at . We present a systematic map of the thermodynamic properties of the LJ spline model from molecular dynamics and Gibbs ensemble Monte Carlo simulations. Results are presented for gas/liquid, liquid/solid and gas/solid coexistence curves, the Joule-Thomson inversion curve, and several other thermodynamic properties. The critical point for the model is estimated to be and , respectively. The triple point is estimated to and . The coexistence densities, saturation pressure, and supercritical isotherms of the LJ spline model were fairly well represented by the Peng-Robison equation of state. We find that Barker-Henderson perturbation theory works less good for the LJ spline than for the LJ model. The first-order perturbation theory overestimates the critical temperature and pressure by about 10% and 90%, respectively. A second-order perturbation theory is not much better. Our assessment is that the mean compressibility approximation gives a poor representation of the second-order perturbation term. Our main conclusion is that we at the moment do not have a theory or model that adequately represents the thermodynamic properties of the LJ spline system.

CO2 sorption in triethyl(butyl)phosphonium 2-cyanopyrrolide ionic liquid via first principles simulations

Predicting gas-liquid solubility using traditional force field based approaches is not reliable in cases where gas molecules can chemically complex with the liquid phase. Using electronic structure calculations to calculate system energy on-the-fly can overcome this deficiency albeit at a higher computational cost. Here, we use first-principles simulations to predict the solubility of CO2 in triethyl(butyl)phosphonium 2-cyanopyrrolide ionic liquid (IL). We utilize Gibbs ensemble Monte Carlo (GEMC) simulations to predict solubility, and employ ab-initio molecular dynamics (AIMD) to study dynamics and structural properties. Solubility prediction from GEMC simulations is in good agreement with the experimental data. Analyzing complexation states of CO2 molecules reveals that approximately 30% of CO2 molecules are chemically bonded with the anion component of the IL. The combined distribution function analysis shows that upon reaction, the OCO angle of the CO2-anion complex varies from 132 to 134 °. Despite these structural changes, plots of mean square displacement (MSD) show minimal difference between the dynamics of neat and IL with dissolved CO2. The vibrational spectra obtained from AIMD simulations match reasonably well for pure and mixed IL as compared to the experimental data.

Separation of nitrogen dioxide from the gas mixture with nitrogen by use of ZIF materials; computer simulation studies

The separation of nitrogen dioxide (NO2) from the gas mixture with nitrogen (N2) by selective adsorption and diffusion in three Zeolitic Imidazolate Frameworks (ZIFs) is examined by Gibbs Ensemble Monte Carlo (GEMC) and Molecular Dynamics (MD) computer simulations at 373 K, which is a temperature that can be found in the exhaust gases of combustion engines and furnaces. This temperature ensures that the chemical equilibrium reaction between NO2 and N2O4 can be neglected since at this temperature almost only NO2 exists. At 373 K, the membrane selectivity of NO2 over N2 is predicted to be about 4.2 (ZIF-90), 5.2 (ZIF-8), and 10.6 (ZIF-78). In addition, it is shown that N2O4 can be neither adsorbed nor diffuses inside these ZIF materials, which makes these ZIFs a promising materials for separation by membranes. Due to molecular sieving, only N2 and NO2 can pass the ZIF membranes.

Cage Occupancies in Nitrogen Clathrate Hydrates from Monte Carlo Simulations

Comparisons of Gibbs ensemble Monte Carlo simulations with experimental data for the cage occupancies in N2 clathrate hydrates are performed to assess the accuracy of such simulations, to refine the effective potentials employed, and to help interpret recently measured large cage over small cage occupancy ratios. Different sets of interaction potentials for N2–N2, N2–H2O, and H2O–H2O interactions are considered. Some of them fail to reproduce the known experimental fact that some large cages are doubly occupied at 273 K and high pressures. The best agreement between simulations and experiments is obtained when using a new N–O interaction potential derived in this work by averaging an ab initio potential energy surface for the N2–H2O dimer.

Molecular simulation of the constant composition expansion experiment in shale multi-scale systems

The phase behavior of hydrocarbons in shale reservoirs has garnered increasing attention in the petroleum industry. Significant differences in the phase behavior of petroleum fluids between conventional reservoirs and shale reservoirs have been observed. Because of the existence of nano-scale porous media in shale reservoirs, there are substantial surface fluid interactions that can lead to a heterogeneous distribution of molecules and an alteration of the fluid phase behavior. In this work, for the first time, we use molecular simulations to investigate the confinement effect on the phase behavior of reservoir fluids in a multi-scale (nano-scale + macro-scale) media. The model has a confined region that mimics the nano-scale porous media and a bulk region that represents the macro-scale porous media (macropores and fractures) in shale reservoirs. For the first time, we use the Gibbs ensemble Monte Carlo (GEMC) simulation at imposed pressures to simulate the constant composition expansion (CCE) experiment for multi-component hydrocarbon mixtures in a multi-scale porous media. Our results show that because of the confinement effect, the composition of the bulk and the confined regions differ significantly. In our molecular simulations, the composition of the lighter components in the bulk region increases, while heavier components are mainly trapped in the confined region. These changes in the composition cause a significant shift or disappearance of the saturation pressure in the bulk region of the multi-component multi-scale hydrocarbon systems.

Cage occupancies, lattice constants, and guest chemical potentials for structure II hydrogen clathrate hydrate from Gibbs ensemble Monte Carlo simulations.

In this paper, equilibrium properties of structure II hydrates of hydrogen were determined from Monte Carlo simulations in the isothermal-isobaric Gibbs ensemble. Water and hydrogen molecules are described by the TIP4P/Ice and Silvera-Goldman models, respectively. The use of the Gibbs ensemble has many key advantages for the simulation of hydrates. By the separation of hydrogen vapor and hydrate phases into their own domains, coupled with transfer moves of hydrogen molecules between domains, cage occupancies were determined. Furthermore, the choice of this ensemble also allows equilibrium lattice constants and guest molecule chemical potentials to be straightforwardly estimated. Results for hydrogen mass fractions indicate reasonable agreement with prior simulation data and theoretical models, while detailed analysis of cage occupancy distributions and neighboring cage pair occupancy combinations gives valuable insight into the behavior of this hydrate at the inter-cage scale. These results will aid in the construction of theoretical models, for which knowledge of the occupancy of neighboring cages is of great importance. In support of previous experimental and theoretical works, we also find evidence of double occupancy of a few small cages inside of the hydrate stability zone, albeit at very high pressures; approximately 0.1% of small cages are doubly occupied at 300 MPa, for temperatures of 225 K and 250 K.

Pressure and Temperature Tuning of Gas-Expanded Liquid Structure and Dynamics.

Carbon dioxide-expanded liquids (CXLs) represent an important class of reaction media that provide tunability of mass transport, solvation, and solubility. Their properties have been demonstrated to provide advantages over traditional organic solvents. However, the molecular-level effects of the CO2 expansion on the structure and dynamics of the liquid that lead to this result have not been fully explored. To address this question, we have used molecular simulations to examine the behavior of two CXLs relevant to the hydroformylation of 1-octene, which has been demonstrated to benefit from the use of gas-expanded reaction media. Specifically, the phase equilibrium properties of CO2-expanded 1-octene and nonanal are calculated as functions of temperature and pressure using Gibbs ensemble Monte Carlo simulations to determine the pressure-composition phase diagrams and volume expansion. In addition, molecular dynamics (MD) simulations were conducted to compute the liquid structure, diffusion coefficients, and shear viscosities. The simulated phase diagrams are in excellent agreement with previous experimental data when available, validating the models used. The MD simulations reveal a direct, linear relationship between the liquid viscosity and the volume expansion, which has not been previously reported. In contrast, deviations from such a relationship are observed for the diffusion coefficient at large volume expansion, indicating that a single Stokes-Einstein relation cannot describe the behavior at all pressures.

Deep neural network learning of complex binary sorption equilibria from molecular simulation data.

We employed deep neural networks (NNs) as an efficient and intelligent surrogate of molecular simulations for complex sorption equilibria using probabilistic modeling. Canonical (N 1 N 2 VT) Gibbs ensemble Monte Carlo simulations were performed to model a single-stage equilibrium desorptive drying process for (1,4-butanediol or 1,5-pentanediol)/water and 1,5-pentanediol/ethanol from all-silica MFI zeolite and 1,5-pentanediol/water from all-silica LTA zeolite. A multi-task deep NN was trained on the simulation data to predict equilibrium loadings as a function of thermodynamic state variables. The NN accurately reproduces simulation results and is able to obtain a continuous isotherm function. Its predictions can be therefore utilized to facilitate optimization of desorption conditions, which requires a laborious iterative search if undertaken by simulation alone. Furthermore, it learns information about the binary sorption equilibria as hidden layer representations. This allows for application of transfer learning with limited data by fine-tuning a pretrained NN for a different alkanediol/solvent/zeolite system.

Prediction of phase equilibria and Gibbs free energies of transfer using molecular exchange Monte Carlo in the Gibbs ensemble

The molecular exchange Monte Carlo (MEMC) method is extended to Gibbs ensemble Monte Carlo (GEMC). The utility of the MEMC method is demonstrated through the calculation of the free energies of transfer of n-alkanes from vapor into liquid 1-octanol, n-hexadecane, and 2,2,4-trimethylpentane using isobaric-isothermal GEMC simulations. Calculations are performed for both the Transferable Potentials for Phase Equilibria (TraPPE) and Mie potentials. For all solutes, the Mie potentials predict free energies of transfer that are within 0.1 kcal/mol of experiment for solvation in n-hexadecane and 2,2,4-trimethylpentane. For TraPPE, solutes shorter than n-butane show similar agreement with experiment, but free energies for n-octane in n-hexadencane and 2,2,4-trimethylpentane are over-predicted by approximately 0.25 and 0.5 kcal/mol, respectively.

Vapor–liquid equilibrium and molecular simulation data for carbon dioxide (CO2) + trans-1,3,3,3-tetrafluoroprop-1-ene (R-1234ze(E)) mixture at temperatures from 283.32 to 353.02 K and pressures up to 7.6 MPa

Isothermal Vapor–Liquid Equilibrium (VLE) data for the binary mixture of CO2 + R-1234ze(E)(trans-1,3,3,3-tetrafluoroprop-1-ene) were measured using a static-analytic method apparatus at seven temperatures between 283.32 and 353.02 K and pressures up to 7.6 MPa. For temperatures over the critical temperature of pure CO2, the critical compositions and pressures of binary mixtures were approximated applying power laws with asymptotic behavior at critical point. The data were well correlated using the Peng–Robinson equation of state incorporating the generalized alpha function, with the Wong–Sandler mixing rules involving NRTL activity coefficient model. The experimental and correlated phase compositions were compared with Gibbs Ensemble Monte Carlo simulation data obtained at 293.15 K and 353.02 K. Though the molecular simulation predictions at low temperature (293.15 K) properly match experimental data, a slight shift appeared at high temperature (353.02 K), and the expected convergence of phase compositions near critical point was not fully observed.

The Influence of Bond Angle on Thermophysical Properties of Three-Center Lennard-Jones Fluids: Computer Simulation and Theory

Abstract We carry out extensive computer simulations to study the phase equilibrium, thermodynamics, and diffusion coefficient of three-center Lennard-Jones (3CLJ) fluids with an emphasis on the effects of bond angle on these properties. We take into account several bond angles ranging from θ 0 = 60 to 180 degrees and two equilibrium bond elongations le = 1.0 and 0.5 (in Lennard-Jones length unit). Moreover, we study the fully flexible (FF) 3CLJ fluids for molecules with bond length 1.0. Gibbs ensemble Monte Carlo (MC) simulations are performed to compute the densities of the vapor-liquid coexisting phases and the vapor pressure, and direct three-phase (vapor-liquid-vapor) molecular dynamics (MD) simulations are carried out to calculate the surface tension. We then apply constant NVT MC simulations to obtain the internal energy, the pressure, and the pair correlation function, and utilize equilibrium MD simulations to compute the diffusion coefficient of systems with le = 1.0. In addition to MD simulations, the modified Cohen-Turnbull (mCT) theory is used to to compute the diffusion coefficient and the mean free volume appearing in the mCT relation is provided by the results of the Generic van der Waals (GvdW) theory. We show that the mCT theory is capable of reproducing the MD simulation values quite well over a wide range of density with slight overestimations at medium range. The angle dependence of different thermophysical properties are analyzed and discussed in details.

Using molecular simulations to probe pore structures and polymer partitioning in size exclusion chromatography

Molecular simulations have been extensively utilized to understand and predict the polymer partitioning in size-exclusion chromatography (SEC). However, idealized pore models (e.g., cylindrical, spherical, and slit pores) were often used to represent the porous media in an SEC column, which leads to significant deviations in describing the geometry and the size of the pores. In this work, several complex pore models were derived from body-centered cubic, random, and gel packing of monodisperse spherical sol particles using simulation methodology. The mechanical stabilities of these structures were determined based on particle coordination numbers. Pore size distributions of these porous structures were compared to a commercially available, wide-pore superficially porous particle. Then, Gibbs ensemble Monte Carlo simulations were performed to compute the pore-to-bulk partitioning coefficient KSEC of a polymer chain with complex pore models. The effects of particle size, packing structure, and porosity on KSEC were explored. In addition, structural analysis provides insight into the conformation of polymers in the pores and its effect on the partitioning behavior. This study promotes the understanding of pore structures in SEC columns and enables more accurate predictions of KSEC with less ambiguity in pore geometry.

Investigating adsorption- and diffusion selectivity of CO2 and CH4 from air on zeolitic imidazolate Framework-78 using molecular simulations

Abstract The adsorption and diffusion of N2, O2 representing air, in ternary mixtures with CO2 or CH4, on Zeolitic Imidazolate Framework-78 are investigated by Gibbs ensemble Monte Carlo (GEMC) simulations for adsorption and by molecular dynamics (MD) simulations for diffusion. The adsorption and diffusion data are used for determining the membrane selectivities. N2/O2 ratios are found to be nearly the same in the gas and the adsorbed phase. At 228 K and low concentrations the CO2/N2 adsorption selectivity is found to attain values of as large as 140, further increasing with decreasing CO2 concentration.