Gibbs Ensemble Monte Carlo Simulations Research Articles

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.

United atom forcefield for vapor-liquid equilibrium (VLE) properties of cyclic and polycyclic compounds from Monte Carlo simulations

Gibbs ensemble Monte Carlo and NPT simulations are applied to the prediction of VLE properties for cyclic and polycyclic molecules with internal flexibility. An extension of the TraPPE-UA forcefield is defined, considering the molecules fully flexible with appropriate intramolecular potential terms (bond stretching, bending, torsion, Van der Waals) in cycles as well as in side chains. New United Atom types are proposed to handle the ternary and quaternary carbons involved in naphthenes and aromatics bearing alkyl groups or naphthenoaromatics, and published TraPPE-UA parameters are kept unchanged for other groups. The extended forcefield is tested for accuracy and internal consistency on a set of 34 compounds comprising between 5 and 28 carbon atoms including naphthenic hydrocarbons, aromatic hydrocarbons, naphthenoaromatics, and thiophenic compounds. This set of compounds contains several multi-functional polycyclic molecules which can now be described with the forcefield extensions provided by this work (e.g. 5α cholestane, 1-hexyl tetralin) and that are representative of a much wider range of compounds for which predictions through theoretical methods are difficult to be achieved. Comparison with available experimental data on normal boiling temperature, vapor pressure, and liquid density, reveals average absolute deviations of 6.4 K on normal boiling temperature and 0.91% on standard liquid density, without noticeable increase of deviations with molecular weight. Computational efficiency allows simulations for compounds containing up to 30 carbon atoms.

Probing Additive Loading in the Lamellar Phase of a Nonionic Surfactant: Gibbs Ensemble Monte Carlo Simulations Using the SDK Force Field.

Understanding solute uptake into soft microstructured materials, such as bilayers and worm-like and spherical micelles, is of interest in the pharmaceutical, agricultural, and personal care industries. To obtain molecular-level insight on the effects of solutes loading into a lamellar phase, we utilize the Shinoda-Devane-Klein (SDK) coarse-grained force field in conjunction with configurational-bias Monte Carlo simulations in the osmotic Gibbs ensemble. The lamellar phase is comprised of a bilayer formed by triethylene glycol mono- n-decyl ether (C10E3) surfactants surrounded by water with a 50:50 surfactant/water weight ratio. We study both the unary adsorption isotherm and the effects on bilayer structure and stability caused by n-nonane, 1-hexanol, and ethyl butyrate at several different reduced reservoir pressures. The nonpolar n-nonane molecules load near the center of the bilayer. In contrast, the polar 1-hexanol and ethyl butyrate molecules both load with their polar bead close to the surfactant head groups. Near the center of the bilayer, none of the solute molecules exhibits a significant orientational preference. Solute molecules adsorbed near the polar groups of the surfactant chains show a preference for orientations perpendicular to the interface, and this alignment with the long axis of the surfactant molecules is most pronounced for 1-hexanol. Loading of n-nonane leads to an increase of the bilayer thickness, but does not affect the surface area per surfactant. Loading of polar additives leads to both lateral and transverse swelling. The reduced Henry's law constants of adsorption (expressed as a molar ratio of additive to surfactant per reduced pressure) are 0.23, 1.4, and 14 for n-nonane, 1-hexanol, and ethyl butyrate, respectively, and it appears that the SDK force field significantly overestimates the ethyl butyrate-surfactant interactions.

Simulation of osmotic pressure across an amorphous semipermeable membrane

Molecular simulations were employed to investigate osmosis for systems far from ideal conditions. Osmotic pressure was calculated for of a fluid mixture separated from a pure fluid by a fictitious semipermeable membrane using Gibbs ensemble Monte Carlo (GEMC) and compared with that of the two phases separated by a physical amorphous polyamide membrane using molecular dynamics. The calculated osmotic pressure for these two simulation set-ups was compared with various predictive models for a range of solute concentrations and solution densities. Both equilibrium-based theories, which are mostly limited in their applicability to near-ideal conditions or low solute concentrations, and non-equilibrium theories that extend to more concentrated solution conditions fared well against the GEMC simulation in their expected range, but did not adequately predict the osmotic pressure calculated according to simulations across the physical membrane for finite system size and under the limiting conditions considered. The discrepancies were particularly large for high solute concentrations and low solution densities. The observed differences are discussed in terms of solution inhomogeneities resulting from particle-membrane interactions.

Understanding the Molecular Weight Dependence of χ and the Effect of Dispersity on Polymer Blend Phase Diagrams

Gibbs ensemble Monte Carlo simulations and cloud point measurements were performed to understand the molecular weight dependence of χ and the effect of dispersity on the phase behavior of polymer mixtures. Oligomeric blends consisting of poly(ethylene-alt-propylene) (PEP) and poly(ethylene oxide) dimethyl ether (PEO) were used as the model systems. First, the molecular weight dependence of χ for PEP/PEO mixtures was studied using simulations and experiments for PEP/PEO mixtures with various molecular weights. An empirical model with a single adjustable parameter kij is used to quantify this molecular weight dependence, and it allows for the accurate prediction of χ of PEP/PEO mixtures with arbitrary molecular weights. Second, the effects of molecular weight distribution (MWD) and dispersity (Đ) of PEO on the PEP/PEO phase diagram were investigated via both simulations and experiments. When PEO is relatively monodisperse (Đ < 1.2), the phase diagram is found to be insensitive to either MWD or Đ, despite di...

A Novel Microemulsion Phase Transition: Toward the Elucidation of Third-Phase Formation in Spent Nuclear Fuel Reprocessing.

We present evidence that the transition between organic and third phases, which can be observed in the plutonium uranium reduction extraction (PUREX) process at high metal loading, is an unusual transition between two isotropic bicontinuous microemulsion phases. As this system contains so many components, however, we have been seeking first to investigate the properties of a simpler system, namely, the related metal-free, quaternary water/n-dodecane/nitric acid/tributyl phosphate (TBP) system. This quaternary system has been shown to exhibit, under appropriate conditions, three coexisting phases: a light organic phase, an aqueous phase, and the so-called third phase. In the current work, we focused on the coexistence of the light organic phase with the third phase. Using Gibbs ensemble Monte Carlo (GEMC) simulations, we found coexistence of a phase rich in nitric acid and dilute in n-dodecane (the third phase) with a phase more dilute in nitric acid but rich in n-dodecane (the light organic phase). The compositions and densities of these two coexisting phases determined using the simulations were in good agreement with those determined experimentally. Because such systems are generally dense and the molecules involved are not simple, the particle exchange rate in their GEMC simulations can be rather low. To test whether a system having a composition between those of the observed third and organic phases is indeed unstable with respect to phase separation, we used the Bennett acceptance ratio method to calculate the Gibbs energies of the homogeneous phase and the weighted average of the two coexisting phases, where the compositions of these phases were taken both from experimental results and from the results of the GEMC simulations. Both demixed states were determined to have statistically significant lower Gibbs energies than the uniform, mixed phase, providing confirmation that the GEMC simulations correctly predicted the phase separation. Snapshots from the simulations and a cluster analysis of the organic and third phases revealed structures akin to bicontinuous microemulsion phases, with the polar species residing within a mesh and with the surface of the mesh formed by amphiphilic TBP molecules. The nonpolar n-dodecane molecules were observed in these snapshots to be outside this mesh. The only large-scale structural differences observed between the two phases were the dimensions of the mesh. Evidence for the correctness of these structures was provided by the results of small-angle X-ray scattering (SAXS) studies, where the profiles obtained for both the organic and third phases agreed well with those calculated from simulations. Finally, we looked at the microscopic structures of the two phases. In the organic phase, the basic motif was observed to be one nitric acid molecule hydrogen-bonded to a TBP molecule. In the third phase, the most common structure was that of the hydrogen-bonded TBP-HNO3-HNO3 chain. A cluster analysis provided evidence for TBP forming an extended, connected network in both phases. Studies of the effects of metal ions on these systems will be presented elsewhere.

Solubility of CO2 in triglycerides using Monte Carlo simulations

The knowledge of solubility of different gases in triglycerides is required to design efficient cell disruption processes for biodiesel production from wet biomass. In this combined experimental and simulation study, we determine carbon dioxide solubility in long (C18) and short (C4) chain triglycerides, and provide microscopic picture of CO2–triglyceride interaction. The solubility of CO2 in tributyrin (a small chain triglyceride) was determined using a pressure drop gas apparatus at P = 800–1700 kPa and T = 283.2–303.2 K. Gibbs ensemble Monte Carlo simulations with configurational biased Monte Carlo algorithm were used to compute the solubility of CO2 in both a long chain triglyceride (to mimic canola oil) and tributyrin at different pressures and temperatures. The transferable potentials for phase equilibria force field was used for modeling triglycerides and CO2. Predicted densities and viscosities of both canola oil and tributyrin were in good agreement with the experimental data. We find that the computed solubilities are in excellent agreement with the experimental results for both CO2–canola oil and CO2–tributyrin where the mean deviation were 2.85% and 9.51%, respectively. The structural analysis of liquid phase indicates that CO2 has slightly higher preference to bind with carbonyl group than the nonpolar segments.

Liquid-Vapor Phase Diagram of RPBE-D3 Water: Electronic Properties along the Coexistence Curve and in the Supercritical Phase.

On the basis of ab initio Gibbs ensemble Monte Carlo simulations, we map the liquid-vapor phase diagram of water described by the RPBE density functional supplemented by D3 dispersion corrections and estimate the critical point by density extrapolation. Knowing the approximate location of the critical point, two sets of ab initio molecular dynamics simulations at gas-like and liquid-like densities deep in the supercritical phase of water are carried out where particular attention is payed to ergodic sampling in view of large correlation lengths and long correlation times. Structural, H-bonding, and dipolar properties of RPBE-D3 water are analyzed along the liquid-vapor coexistence curve upon moving toward the critical point and compared to those in the supercritical state. The properties of high-density supercritical water are astonishingly similar to those of the liquid on the coexistence curve under subcritical conditions at comparable density. Upon decomposing the molecular dipole moments into purely configurational and electronic polarization/charge-transfer contributions, it is demonstrated that the latter play a decreasing role in liquid water upon approaching the critical point on the coexistence curve. Moreover, these many-body effects are systematically suppressed in supercritical water due to the significantly reduced H-bonding network.

Virial coefficients, equation of state, and demixing of binary asymmetric nonadditive hard-disk mixtures

Values of the fifth virial coefficient, compressibility factors, and fluid-fluid coexistence curves of binary asymmetric nonadditive mixtures of hard disks are reported. The former correspond to a wide range of size ratios and positive nonadditivities and have been obtained through a standard Monte Carlo method for the computation of the corresponding cluster integrals. The compressibility factors as functions of density, derived from canonical Monte Carlo simulations, have been obtained for two values of the size ratio (q = 0.4 and q = 0.5), a value of the nonadditivity parameter (Δ = 0.3), and five values of the mole fraction of the species with the biggest diameter (x1 = 0.1, 0.3, 0.5, 0.7, and 0.9). Some points of the coexistence line relative to the fluid-fluid phase transition for the same values of the size ratios and nonadditivity parameter have been obtained from Gibbs ensemble Monte Carlo simulations. A comparison is made between the numerical results and those that follow from some theoretical equations of state.

NLDFT Pore Size Distribution in Amorphous Microporous Materials.

The pore size distribution (PSD) is one of the most important properties when characterizing and designing materials for gas storage and separation applications. Experimentally, one of the current standards for determining microscopic PSD is using indirect molecular adsorption methods such as nonlocal density functional theory (NLDFT) and N2 isotherms at 77 K. Because determining the PSD from NLDFT is an indirect method, the validation can be a nontrivial task for amorphous microporous materials. This is especially crucial since this method is known to produce artifacts. In this work, the accuracy of NLDFT PSD was compared against the exact geometric PSD for 11 different simulated amorphous microporous materials. The geometric surface area and micropore volumes of these materials were between 5 and 1698 m2/g and 0.039 and 0.55 cm3/g, respectively. N2 isotherms at 77 K were constructed using Gibbs ensemble Monte Carlo (GEMC) simulations. Our results show that the discrepancies between NLDFT and geometric PSD are significant. NLDFT PSD produced several artificial gaps and peaks that were further confirmed by the coordinates of inserted particles of a specific size. We found that dominant peaks from NLDFT typically reported in the literature do not necessarily represent the truly dominant pore size within the system. The confirmation provides concrete evidence for artifacts that arise from the NLDFT method. Furthermore, a sensitivity analysis was performed to show the high dependency of PSD as a function of the regularization parameter, λ. A higher value of λ produced a broader and smoother PSD that closely resembles geometric PSD. As an alternative, a new criterion for choosing λ, called here the smooth-shift method (SSNLDFT), is proposed that tuned the NLDFT PSD to better match the true geometric PSD. Using the geometric pore size distribution as our reference, the smooth-shift method reduced the root-mean-square deviation by ∼70% when the geometric surface area of the material is greater than 100 m2/g.

Determining the critical points and coexistence curve of multi-species equilibria using Gaussian Process Regression

We present a method for reconstructing the liquid-vapor coexistence curve and associated location of critical points of multi-species mixtures from experimental or simulation data, using Gaussian Process Regression. The distinguishing feature of the proposed method is that it does not require knowledge of the equation of state, as assumed by the currently prevalent methods for calculating the coexistence properties of mixtures. The proposed technique is validated using synthetic data generated from the Peng-Robinson equation of state, experimental data, as well as Gibbs Ensemble Monte Carlo simulation data. These validation studies show that the proposed method is robust to statistical noise and scarcity of data even in the vicinity of the critical point.

Molecular Simulation Study of the Performance of Supported Ionic Liquid Phase Materials for the Separation of Carbon Dioxide from Methane and Hydrogen

Molecular dynamics and Gibbs ensemble Monte Carlo simulations were used to compute the self-diffusion coefficients and solubilities of CO2, CH4, and H2 in model membranes consisting of slit pores with diameters of 2 and 5 nm. Solubility selectivities, diffusion selectivities, and permselectivities of CO2 for binary gas mixtures of CO2/CH4 and CO2/H2 were also computed. The calculations were repeated for the same pores filled with the ionic liquid (IL) 1-n-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([C4mim]+[Tf2N]−) and for bulk IL. The bulk IL system was used as a model for a supported ionic liquid membrane separator having large pores, while the confined IL systems were used to assess whether extreme nanoconfinement of ILs has an effect on permselectivity. Permselectivities were about a factor of 10 higher in all the IL systems compared to the empty nanopores. Nanoconfinement tends to increase the solubility and decrease the diffusivity of all gases relative to the bulk IL. The bulk IL h...

Bis-chlorinated aromatics adsorption in Faujasites investigated by molecular simulation-influence of Na+ cation

Molecular simulations have been employed in order to explore at the microscopic scale the adsorption of bis-chlorinated aromatics (1,2- and 1,3-dichlorobenzene) in Faujasite, a large pore openings zeolite. Both, the purely siliceous and the sodium cation exchanged forms of zeolite structures have been respectively considered, aiming to clearly determine the role of the charge compensating cation in the adsorption phenomenon of the investigated aromatics. A suited combination of classical and electronic structure simulation tools provided a clear overall picture of the adsorption process, from both local and global points of view, well matching with the accessible experimental data. The adsorbate preferential locations, adsorbate/zeolite interaction nature and geometry as well as the adsorption energy were extracted from Density Functional Theory calculations. Furthermore, on the basis of classical force fields, the Gibbs ensemble Monte Carlo simulations allowed predicting the room temperature (298 K) adsorption isotherms for the investigated molecules in Faujasite, within the purely siliceous and Na+ exchanged form, with mobile extra-framework cations upon the adsorption process. Finally, we accomplished a detailed analysis of the microscopic mechanism in play along the whole adsorption process, with a special highlight to the understanding of the interaction geometry of the molecule with the sodium cation, in function of its crystallographic site. Location of the charge compensating cation has been found to influence the nature of interaction with the adsorbate molecule.

Uncertainty quantification and propagation of errors of the Lennard-Jones 12-6 parameters for n-alkanes.

Molecular simulation has the ability to predict various physical properties that are difficult to obtain experimentally. For example, we implement molecular simulation to predict the critical constants (i.e., critical temperature, critical density, critical pressure, and critical compressibility factor) for large n-alkanes that thermally decompose experimentally (as large as C48). Historically, molecular simulation has been viewed as a tool that is limited to providing qualitative insight. One key reason for this perceived weakness in molecular simulation is the difficulty to quantify the uncertainty in the results. This is because molecular simulations have many sources of uncertainty that propagate and are difficult to quantify. We investigate one of the most important sources of uncertainty, namely, the intermolecular force field parameters. Specifically, we quantify the uncertainty in the Lennard-Jones (LJ) 12-6 parameters for the CH4, CH3, and CH2 united-atom interaction sites. We then demonstrate how the uncertainties in the parameters lead to uncertainties in the saturated liquid density and critical constant values obtained from Gibbs Ensemble Monte Carlo simulation. Our results suggest that the uncertainties attributed to the LJ 12-6 parameters are small enough that quantitatively useful estimates of the saturated liquid density and the critical constants can be obtained from molecular simulation.

Importance of ZIF-90 Lattice Flexibility on Diffusion, Permeation, and Lattice Structure for an adsorbed H2/CH4 Gas Mixture: A Re-Examination by Gibbs Ensemble Monte Carlo and Molecular Dynamics Simulations

Adsorption and diffusion of the gas mixture H2/CH4 in the metal–organic framework (MOF) of structure type zeolitic imidazolate framework-90 (ZIF-90) are revisited. While the adsorption can successfully be examined in Gibbs ensemble Monte Carlo (GEMC) simulations using the common approximation of a rigid lattice, the dynamics of methane in ZIF-90 is remarkably influenced by the lattice flexibility. Molecular dynamics (MD) simulations not only show a strong influence of the lattice flexibility on the diffusion of methane but even find a slight structural phase transition of the lattice. This structural change appeared at higher temperatures and was not caused by the content of guest molecules like in most former discoveries of gate opening. For prediction of mixed gas ZIF-90 membrane selectivity, the adsorption and diffusion results show that the high CH4 adsorption selectivity is overcompensated by the high H2 mobility. The comparison of our results for the H2/CH4 membrane selectivity with experimental fin...

Optimised Mie potentials for phase equilibria: application to alkynes

ABSTRACTA transferable united-atom force field, based on Mie potentials, is presented for alkynes. The performance of the optimised Mie potential parameters is assessed for 1-alkynes and 2-alkynes using grand canonical histogram-reweighting Monte Carlo simulations. For each compound, vapour–liquid coexistence curves, vapour pressures, heats of vapourisation, critical properties and normal boiling points are predicted and compared to experiment. Experimental saturated liquid densities are reproduced to within 2% average absolute deviation (AAD), except for 1-hexyne, which are reproduced with 3% AAD. Experimental saturated vapour pressures are reproduced to within 3% AAD, except for 1-pentyne, 2-pentyne and 2-hexyne, where deviations from experiment of up to 20% AAD were observed. Binary phase diagrams, predicted from Gibbs ensemble Monte Carlo simulations, for propane + propyne, propene + propyne and propadiene + propyne, are in close agreement with experiment.

Experimental and Molecular Modeling Study of Bubble Points of Hydrocarbon Mixtures in Nanoporous Media

The shale play resources have played a key role in increasing oil production in the past decade in the United States. The sizes of pores in shales storing the oil are believed to be on the order of nanometers. It is believed that the fluids present in such small nanometer-scale pores have different properties compared to properties measured in the bulk. Fluid saturation pressures at given temperatures, bubble points for oils and dew points for condensates, in the nanopores are affected by the influence of pore walls in the vicinity of the fluid molecules. An approach to bubble point or dew point influences the proportion of liquid or gas produced from a given well and, thus, impacts the economic viability. Hence, an accurate measure of saturation pressures is important. In this paper, we describe experiments in well-characterized synthesized mesoporous materials and present Gibbs ensemble Monte Carlo (GEMC) simulations for understanding the possible reason(s) for observations made in the experiments. The ...

Assessment and Optimization of Configurational-Bias Monte Carlo Particle Swap Strategies for Simulations of Water in the Gibbs Ensemble.

Particle swap moves between phases are usually the rate-limiting step for Gibbs ensemble Monte Carlo (GEMC) simulations of fluid phase equilibria at low reduced temperatures because the acceptance probabilities for these moves can become very low for molecules with articulated architecture and/or highly directional interactions. The configurational-bias Monte Carlo (CBMC) technique can greatly increase the acceptance probabilities, but the efficiency of the CBMC algorithm is influenced by multiple parameters. In this work we assess the performance of different CBMC strategies for GEMC simulations using the SPC/E and TIP4P water models at 283, 343, and 473 K, demonstrate that much higher acceptance probabilities can be achieved than previously reported in the literature, and make recommendations for CBMC strategies leading to optimal efficiency.

Monte Carlo predictions for vapor-liquid equilibria of N-Methyl-2-pyrrolidone + benzene + thiophene

Abstract An extension of the transferable potentials for phase equilibria united-atom (TraPPE-UA) force field to N-Methyl-2-pyrrolidone (NMP) is proposed in this work. The extensional force field can accurately predict the saturated liquid densities, vapor pressures, normal boiling point, heats of vaporization and critical properties of NMP. Gibbs ensemble Monte Carlo (GEMC) simulations are carried out to calculate the isobaric vapor–liquid equilibria for the binary mixtures of benzene + NMP, thiophene + NMP and their ternary mixture at 101.3 kPa. All of these simulation results are in good accordance with the available experimental data without introducing additional binary interaction parameters. In order to study the mechanism of extracting thiophene from benzene by NMP, NPT Monte Carlo simulations are performed to have obtained the radial distribution functions for the ternary liquid mixture benzene + thiophene + NMP. The results show that the oxygen atom of NMP forms stronger hydrogen bond with the hydrogen atom of thiophene than with that of benzene.