The series of metal–organic frameworks M-MOF-74 gained popularity in the field of capture and separation of CO₂ due to the presence of numerous, highly reactive open-metal sites. The description of effective interactions between guest molecules and open-metal sites without accounting for polarization effects is challenging but it can significantly reduce the computational cost of simulations. In this study, we propose a non-polarizable force field for CO₂, and H₂ adsorption in M-MOF-74 (M = Ni, Cu, Co, Fe, Mn, Zn) by scaling the Coulombic interactions of M-MOF-74 atoms, and Lennard-Jones interaction potentials between the center of mass of H₂ and the open-metal centers. The presented force field is based on UFF and DREIDING parameters, characterized by high transferability and efficiency. The quantum behavior of H₂ at cryogenic temperatures is considered by incorporating Feynman–Hibbs quantum corrections. To validate the force field, the experimental isotherms of CO₂ at 298 K and 10⁻¹ – 10²kPa, the isotherms of H₂ at 77 K and 10⁻⁵ – 10²kPa, the corresponding enthalpy of adsorption, and the binding geometries in the M-MOF-74 series were reproduced using Monte Carlo simulations in the grand-canonical ensemble. The computed loadings, heats of CO₂ and H₂ adsorption, and binding geometries in M-MOF-74 are in very good agreement with the experimental values. The temperature transferability of the force field from 77 K to 87 K, and 298 K was shown for adsorption of H₂. The validated force field was used to study the adsorption and separation of CO₂/H₂ mixtures at 298 K. The adsorption of H₂ practically does not occur when CO₂ is present in the mixture. As indicated from simulated breakthrough curves, the breakthrough time of CO₂ in M-MOF-74 follows the same order as the uptake and the heat of CO₂ adsorption: Ni ¿ Co ¿ Fe ¿ Mn ¿ Zn ¿ Cu. Increasing the feed mole fraction of CO₂ in the breakthrough simulations from 0.1 to 0.9 speeds up the saturation of the adsorbent, leading to a faster exit of CO₂ with the column effluent. The application of the non-polarizable force field allows full investigation of the capture and separation of CO₂ in M-MOF-74, and can be expanded to study multi-component mixtures or industrial reactions in future research.

@en

We study important aspects of shape selectivity effects of zeolites for hydroisomerization of linear alkanes, which produces a myriad of isomers, particularly for long chain hydrocarbons. To investigate the conditions for achieving an optimal yield of branched hydrocarbons, it is important to understand the role of chemical equilibrium in these reversible reactions. We conduct an extensive analysis of shape selectivity effects of different zeolites for the hydroisomerization of C₇ and C₈ isomers at chemical reaction equilibrium conditions. The reaction ensemble Monte Carlo method, coupled with grand-canonical Monte Carlo simulations, is commonly used for computing reaction equilibrium of heterogeneous reactions. The computational demands become prohibitive for a large number of reactions. We used a faster alternative in which reaction equilibrium is obtained by imposing chemical equilibrium in the gas phase and phase equilibrium between the gas phase components and the adsorbed phase counterparts. This effectively mimics the chemical equilibrium distribution in the adsorbed phase. Using Henry’s law at infinite dilution and mixture adsorption isotherm models at elevated pressures, we calculate the adsorbed loadings in the zeolites. This study shows that zeolites with cage or channel-like structures exhibit significant differences in selectivity for alkane isomers. We also observe a minimal impact of pressure on the gas-phase equilibrium of these reactions at typical experimental reaction temperatures 400 − 700 K . This study marks initial strides in understanding the reaction product distribution for long-chain alkanes.

@en

We present the RUPTURA code (https://github.com/iraspa/ruptura) as a free and open-source software package (MIT license) for (1) the simulation of gas adsorption breakthrough curves, (2) mixture prediction using methods like the Ideal Adsorption Solution Theory (IAST), segregated-IAST and explicit isotherm models, and (3) fitting of isotherm models on computed or measured adsorption isotherm data. The combination with the RASPA software enables computation of breakthrough curves directly from adsorption simulations in the grand-canonical ensemble. RUPTURA and RASPA have similar input styles. IAST is implemented near machine precision but we also provide several explicit mixture prediction methods that are non-iterative and potentially faster than IAST. The code supports a wide variety of isotherm models like Langmuir, Anti-Langmuir, BET, Henry, Freundlich, Sips, Langmuir-Freundlich, Redlich-Peterson, Toth, Unilan, O'Brian & Myers, Asymptotic Temkin, and Bingel & Walton. The isotherm model parameters can easily be obtained by the fitting module. Breakthrough plots and animations of the column properties are automatically generated. In addition to highlighting the code, we also review all the developed techniques from literature for mixture prediction, breakthrough simulations, and isotherm model fitting, and provide a tutorial discussing the workflows.

@en

Formic acid production from CO₂ allows the reduction of carbon dioxide emissions while synthesizing a product with a wide range of applications. CO₂ hydrogenation is challenging due to the cost of transition metal catalysts and the toxicity of the transition elements. In this work, the thermodynamic confinement effects of the metal–organic framework UiO-66 on the CO₂ hydrogenation to formic acid were studied by force field-based molecular simulations. The confinement effects of UiO-66 and the metal–organic frameworks Cu-BTC, and IRMOF-1 were compared, to assess the impact of different pore size and metal centers on the production of HCOOH. Monte Carlo simulations in the grand-canonical ensemble were performed in the frameworks, using gas phase mole fractions of CO₂, H₂, and HCOOH at chemical equilibrium, obtained from Continuous Fractional Component Monte Carlo simulations in the Reaction Ensemble. The adsorption isobars of the components in metal–organic frameworks were computed at 298.15 – 800 K, 1 – 60 bar. The enhancement of HCOOH production due to preferential adsorption of HCOOH in metal–organic frameworks was calculated for all studied conditions. UiO-66, Cu-BTC, and IRMOF-1 affect CO₂ hydrogenation reaction, shifting the thermodynamical equilibrium toward HCOOH formation. The prevailing factor is the type of metal center in the metal–organic framework. The confinement effect of Cu-BTC turns out to exceed the enhancement caused by UiO-66, and IRMOF-1. The resulting mole fraction of HCOOH increased by ca. 2000 times compared to the gas phase at 298.15 K, 60 bar. Cu-BTC can be considered as an alternative to improve the production of HCOOH due to elimination of the high-cost temperature elevation, cost reduction of downstream processing methods, and comparable final concentration of HCOOH to the reported concentrations of formate obtained using transition metal catalysts.

@en

Ideal Adsorbed Solution Theory (IAST) is a common method for modelling mixture adsorption isotherms based on pure component isotherms. When the adsorbent has distinct adsorption sites, the segregated version of IAST (SIAST) provides improved adsorbed loadings compared to IAST. We have adopted the concept of SIAST and applied it to an explicit isotherm model which takes into account the different sizes of the adsorbates: the so called Segregated Explicit Isotherm (SEI). The purpose of SEI is to have an explicit adsorption model that can consider both size-effects of the co-adsorbed molecules and surface heterogeneities. In sharp contrast to IAST and SIAST, no iterative scheme is required in case of SEI, which leads to much faster simulations. A comparative study has been performed to analyse the adsorption isotherms calculated using these three methods. The adsorbed loadings predicted by SEI and SIAST are in excellent agreement with the Grand-Canonical Monte Carlo (GCMC) simulation data. The loadings estimated by IAST show considerable deviations from the GCMC data at high pressures. Breakthrough curve modelling is used to compare the effects of these three models at dynamic conditions. The explicit model (SEI) leads to the fastest simulation run time, followed by SIAST.

@en

In this work, adsorption properties of the UiO-66 metal–organic framework were investigated, with particular emphasis on the influence of structural defects. A series of UiO-66 samples were synthesized and characterized using a wide range of experimental techniques. Type I adsorption isotherms for low-temperature adsorption of N₂ and Ar showed that micropore volume and specific surface area significantly increase with the number of defects. Adsorption of hexane isomers in UiO-66 was studied by means of quasi-equilibrated temperature-programmed desorption and adsorption (QE-TPDA) experimental and Monte Carlo simulation techniques. QE-TPDA profiles revealed that only defect-free UiO-66 exhibits distinct two adsorption states. This technique also yielded high-quality adsorption isobars that were successfully recreated using Grand-Canonical Monte Carlo molecular simulations, which, however, required refinement of the existing force fields. The calculations demonstrated the detailed mechanism of adsorption and separation of hexane isomers in the UiO-66 structure. The preferred tetrahedral cages provide suitable voids for bulky molecules, which is the reason for unusual “reverse” selectivity of UiO-66 towards di-branched alkanes. Interconnection of the tetrahedral cavities due to missing organic linkers greatly reduces the selectivity of the defected material.

@en

A reactor model that deconvolutes thermodynamics of adsorption of hydrocarbon in the pores of zeolite Beta, obtained by Configurational-bias Monte Carlo simulations, from intrinsic, intraporous kinetics of hydroisomerization and hydrocracking reactions, provides a good quantitative description of all significant reactions in the kinetic network for interconversion and cracking of different heptane isomers. Activation enthalpies obtained for intraporous reactions follow the expected order according to the carbenium ion formalism: methyl shift< ethyl shift < isom(B) ∼ crack(B2) < crack(B1) < crack(C) ∼ crack(D) < crack(E) and apparently within each isomerization class, in terms of carbenium ions formally involved: sec → tert < sec → sec ∼ tert → tert < tert → sec. except for the ethyl shift reaction forming 3-ethylpentane. Cracking happens primarily through 2,4-dimethylpentane (type B2), regardless of the initial reactant. The model can be subsequently used to separate the effect of pore structure on selective adsorption and on intraporous reaction kinetics. Zeolite Beta will serve as a base case for a comparison of different zeolite structures.

@en

One of the important parameters in water management of proton exchange membranes is the electro-osmotic drag (EOD) coefficient of water. The value of the EOD coefficient is difficult to justify, and available literature data on this for Nafion membranes show scattering from in experiments and simulations. Here, we use a classical all-atom model to compute the EOD coefficient and thermodynamic properties of water from molecular dynamics simulations for temperatures between 330 and 420 K, and for different water contents between λ = 5 and λ = 20. λ is the ratio between the moles of water molecules to the moles of sulfonic acid sites. This classical model does not capture the Grotthuss mechanism; however, it is shown that it can predict the EOD coefficient within the range of experimental values for λ = 5 where the vehicular mechanism dominates proton transfer. For λ > 5, the Grotthuss mechanism becomes dominant. To obtain the EOD coefficient, average velocities of water and ions are computed by imposing different electric fields to the system. Our results show that the velocities of water and hydronium scale linearly with the electric field, resulting in a constant ratio of ca. 0.4 within the error bars. We find that the EOD coefficient of water linearly increases from 2 at λ = 5 to 8 at λ = 20 and the results are not sensitive to temperature. The EOD coefficient at λ = 5 is within the range of experimental values, confirming that the model can capture the vehicular transport of protons well. At λ = 20, due to the absence of proton hopping in the model, the EOD coefficient is overestimated by a factor of 3 compared to experimental values. To analyze the interactions between water and Nafion, the partial molar enthalpies and partial molar volumes of water are computed from molecular dynamics simulations. At different water uptakes, multiple linear regression is used on raw simulation data within a narrow composition range of water inside the Nafion membrane. The partial molar volumes and partial molar excess enthalpies of water asymptotically approach the molar volumes and molar excess enthalpies of pure water for water uptakes above 5. This confirms the model can capture the bulklike behavior of water in the Nafion at high water uptakes.

@en

We present several new major features added to the Monte Carlo (MC) simulation code Brick-CFCMC for phase- and reaction equilibria calculations (https://gitlab.com/ETh_TU_Delft/Brick-CFCMC). The first one is thermodynamic integration for the computation of excess chemical potentials (μex). For this purpose, we implemented the computation of the ensemble average of the derivative of the potential energy with respect to the scaling factor for intermolecular interactions (⟨∂U∂λ⟩). Efficient bookkeeping is implemented so that the quantity ∂U∂λ is updated after every MC trial move with negligible computational cost. We demonstrate the accuracy and reliability of the calculation of μex for sodium chloride in water. Second, we implemented hybrid MC/MD translation and rotation trial moves to increase the efficiency of sampling of the configuration space. In these trial moves, short Molecular Dynamics (MD) trajectories are performed to collectively displace or rotate all molecules in the system. These trajectories are accepted or rejected based on the total energy drift. The efficiency of these trial moves can be tuned by changing the time step and the trajectory length. The new trial moves are demonstrated using MC simulations of a viscous fluid (deep eutectic solvent).

@en

Force field-based molecular simulations were used to calculate thermal expansivities, heat capacities, and Joule-Thomson coefficients of binary (standard) hydrogen-water mixtures for temperatures between 366.15 and 423.15 K and pressures between 50 and 1000 bar. The mole fraction of water in saturated hydrogen-water mixtures in the gas phase ranges from 0.004 to 0.138. The same properties were calculated for pure hydrogen at 323.15 K and pressures between 100 and 1000 bar. Simulations were performed using the TIP3P and a modified TIP4P force field for water and the Marx, Vrabec, Cracknell, Buch, and Hirschfelder force fields for hydrogen. The vapor-liquid equilibria of hydrogen-water mixtures were calculated along the melting line of ice Ih, corresponding to temperatures between 264.21 and 272.4 K, using the TIP3P force field for water and the Marx force field for hydrogen. In this temperature range, the solubilities and the chemical potentials of hydrogen and water were obtained. Based on the computed solubility data of hydrogen in water, the freezing-point depression of water was computed ranging from 264.21 to 272.4 K. The modified TIP4P and Marx force fields were used to improve the solubility calculations of hydrogen-water mixtures reported in our previous study [ Rahbari, A.;et al. J. Chem. Eng. Data 2019, 64, 4103-4115 ] for temperatures between 323 and 423 K and pressures ranging from 100 to 1000 bar. The chemical potentials of ice Ih were calculated as a function of pressure between 100 and 1000 bar, along the melting line for temperatures between 264.21 and 272.4 K, using the IAPWS equation of state for ice Ih. We show that at low pressures, the presence of water has a large effect on the thermodynamic properties of compressed hydrogen. Our conclusions may have consequences for the energetics of a hydrogen refueling station using electrochemical hydrogen compressors.

@en

Vapor-liquid equilibria of xylenes were computed using Monte Carlo simulations in the Gibbs ensemble. For binary mixtures, the predicted composition of the liquid phase is in agreement with experiments. The computed vapor phase densities of each isomer showed an effect on the predicted composition of the vapor phase of the mixture. Unfortunately, the Monte Carlo code used to calculate the vapor-liquid equilibrium contained an error in the acceptance rule for molecule transfer in multi-component mixtures. The Continuous Fractional Component (CFCMC) [1–3] algorithm in the Gibbs ensemble adds an extra molecule -the fractional molecule-per molecule type. In the CFCMC algorithm, trial moves are attempted to transform a fractional molecule of type a in box i to a whole molecule and, simultaneously, transform a molecule of type a in box j [Formula presented] i to a fractional molecule. These trial moves should be accepted according to Ref. [1]: [Formula presented] where [Formula presented] is the number of molecules of type a in box j, [Formula presented] is the number of molecules of type a in box i, and [Formula presented] is the energy difference between the old and the new configuration. Unfortunately, what was implemented (incorrectly) in the computer code is: [Formula presented] where [Formula presented] and [Formula presented] are the total number of molecules of all types in box i and j, respectively, and not only of type a. For the simulation of a single component, [Formula presented] and [Formula presented]. Therefore, the mistake in the acceptance rule only affects the simulations that include more than one component. As such, Figs. 1, Fig. 2, Fig. 3, Fig. 4, and Table 1 of the article are not affected by this error. This error was hard to detect for xylene isomers, as the boiling point of these isomers are close, and therefore [Formula presented] is close to [Formula presented]. The corrected code has been tested to yield the same results as the RASPA software [4,5] for various multi-component systems. The phase composition diagram using the correct code is shown in Fig. 5. The observations based on the incorrect results remain. The composition of the liquid phase is not affected by the choice of force field or method for electrostatic interactions. For all the force fields, the correct vapor phase composition calculated with the Ewald and the Wolf methods are closer than in the incorrect simulations. In the incorrect results, the largest difference between vapor phase composition calculated with the Wolf and the Ewald method was 0.089 [mol/mol]. In the correct results, this difference is reduced to 0.067 [mol/mol]. [Figure presented] The predicted composition of the liquid phase is in excellent agreement with the experimental data. The correct simulations of the vapor phase compositions are closer to the experimental data than the incorrect simulations. The simulations with the incorrect acceptance rules showed an azeotrope behaviour that is not present in the correct simulations of the phase compositions. [Figure presented] The phase composition of the binary mixture and the excess chemical potential are not related. This is consistent with the observations from the incorrect simulations. [Formula presented] The predicted composition of the liquid phase is in excellent agreement with the experimental data. The correct simulations of the vapor phase compositions are closer to the experimental data than the incorrect simulations. The simulations with the incorrect acceptance rules showed an azeotrope behaviour that is not present in the correct simulations. [Figure presented]@en

This article explores how well vapor-liquid equilibria of pure components and binary mixtures of xylenes can be predicted using different force fields in molecular simulations. The accuracy of the Wolf method and the Ewald summation is evaluated. Monte Carlo simulations in the Gibbs ensemble are performed at conditions comparable to experimental data, using four different force fields. Similar results using the Wolf and the Ewald methods can be obtained for the prediction of densities and the phase compositions of binary mixtures. With the Wolf method, up to 50% less CPU time is used compared to the Ewald method, at the cost of accuracy and additional parameter calibration. The densities of p-xylene and m-xylene can be well estimated using the TraPPE-UA and AUA force fields. The largest differences of VLE with experiments are observed for o-xylene. The p-xylene/o-xylene binary mixtures at 6.66 and 81.3 kPa are simulated, leading to an excellent agreement in the predictions of the composition of the liquid phase compared to experiments. The composition of the vapor phase is dominated by the properties of the component with the largest mole fraction in the liquid phase. The accuracy of the predictions of the phase composition are related to the quality of the density predictions of the pure component systems. The phase composition of the binary system of xylenes is very sensitive to slight differences in vapor phase density of each xylene isomer, and how well the differences are captured by the force fields.

@en

Molecular simulations are an excellent tool to study adsorption and diffusion in nanoporous materials. Examples of nanoporous materials are zeolites, carbon nanotubes, clays, metal-organic frameworks (MOFs), covalent organic frameworks (COFs) and zeolitic imidazolate frameworks (ZIFs). The molecular confinement these materials offer has been exploited in adsorption and catalysis for almost 50 years. Molecular simulations have provided understanding of the underlying shape selectivity, and adsorption and diffusion effects. Much of the reliability of the modeling predictions depends on the accuracy and transferability of the force field. However, flexibility and the chemical and structural diversity of MOFs add significant challenges for engineering force fields that are able to reproduce experimentally observed structural and dynamic properties. Recent developments in design, parameterization, and implementation of force fields for MOFs and zeolites are reviewed.

@en

Computational chemistry uses computer simulation to assist in solving chemical problems. Typical workflows of computational chemists include the use of dozens of utilities. 3D modeling programs are powerful tools that help researchers visualize their work and create illustrative graphics. In this review, we describe and highlight tools and visualization packages that are commonly used in the field of (bio-)chemistry and material science.

@en

The separation of light olefins from paraffins via cryogenic distillation is a very energy intensive process. Solid adsorbents and especially metal-organic frameworks with open metal sites have the potential to significantly lower the required energy. Specifically, M-MOF-74 has drawn considerable attention for application in olefin/paraffin separation. To investigate how the separation proceeds on a molecular level and to design better materials, molecular simulation can be a useful tool. Unfortunately, it is still a challenge to model the adsorption behavior of many adsorbates in metal-organic frameworks with open metal sites. Previously, the inclusion of explicit polarization has been suggested to improve the quality of classical force fields for such systems. Here, the potential of polarizable force fields for the description of olefins and paraffins in metal-organic frameworks with open metal sites is investigated. In particular, heats of adsorption, binding geometries, and adsorption isotherms are calculated for C2H4, C2H6, C3H6, and C3H8 in M-MOF-74 (with M = Co, Mn, Fe, and Ni). In this study, no force field parameters are adjusted to improve the model. The results show that including explicit polarization significantly improves the description of the adsorption in comparison to non-polarizable generic force fields which do not consider explicit polarization. The study also reveals that simulation predictions are sensitive to the assigned repulsive potential and framework charges. A fully re-parametrized polarizable force field may have the capability to improve the predictions even further.

@en

Many important industrial separation processes based on adsorption operate close to saturation. In this regime, the underlying adsorption processes are mostly driven by entropic forces. At equilibrium, the entropy of adsorption is closely related to the enthalpy of adsorption. Thus, studying the behavior of the enthalpy of adsorption as a function of loading is fundamental in the understanding of the separation processes. Unfortunately, close to saturation, the enthalpy of adsorption is hard to measure experimentally and hard to compute in simulations. In simulations the enthalpy of adsorption is usually obtained from energy/particle fluctuations in the grand-canonical ensemble, but this methodology is hampered by vanishing insertions/deletions at high loading. To investigate the fundamental behavior of the enthalpy and entropy of adsorption at high loading, we develop a simplistic model of adsorption in a channel and show that at saturation the enthalpy of adsorption diverges to large positive values due to repulsive intermolecular interactions. However, there are many systems that can avoid repulsive intermolecular interactions and hence do not show this drastic increase in enthalpy of adsorption close to saturation. We find that the conventional grand-canonical Monte Carlo (GCMC) method is incapable of determining the enthalpy of adsorption from energy/particle fluctuations at high loading. Here, we show that by using the Continuous Fractional Component Monte Carlo (CFCMC), the enthalpy of adsorption close to saturation conditions can be reliably obtained from the energy/particle fluctuations in the grand-canonical ensemble. The best method to study properties at saturation is the NVT energy (local-) slope methodology. @en

Monte Carlo simulations are used to calculate the solubility of natural gas components in ionic liquids (ILs) and Selexol, which is a mixture of poly(ethylene glycol) dimethyl ethers. The solubility of the pure gases carbon dioxide (CO₂), methane (CH₄), ethane (C₂H₆), and sulfur dioxide (SO₂) in the ILs 1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([C_nmim][Tf₂N], n = 4, 6), 1-ethyl-3-methylimidazolium diethylphosphate ([emim][dep]), and Selexol (CH₃O[CH₂CH₂O]_nCH₃, n = 4, 6) have been computed at 313.15 K and several pressures. The gas solubility trend observed in the experiments and simulations is: SO₂ > CO₂ > C₂H₆ > CH₄. Overall, the Monte Carlo simulation results are in quantitative agreement with existing experimental data. Molecular simulation is an excellent tool to predict gas solubilities in solvents and may be used as a screening tool to navigate through the large number of theoretically possible ILs.

@en