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).

We introduce an alternative method to perform free energy calculations for mixtures at multiple temperatures and pressures from a single simulation, by combining umbrella sampling and the continuous fractional component Monte Carlo method. One can perform a simulation of a mixture at a certain pressure and temperature and accurately compute the chemical potential at other pressures and temperatures close to the simulation conditions. This method has the following advantages: (1) Accurate estimates of the chemical potential as a function of pressure and temperature are obtained from a single state simulation without additional postprocessing. This can potentially reduce the number of simulations of a system for free energy calculations for a specific temperature and/or pressure range. (2) Partial molar volumes and enthalpies are obtained directly from the estimated chemical potentials. We tested our method for a Lennard-Jones system, aqueous mixtures of methanol at T = 298 K and P = 1 bar, and a mixture of ammonia, nitrogen, and hydrogen at T = 573 K and P = 800 bar. For pure methanol (N = 410 molecules), we observed that the estimated chemical potentials from umbrella sampling are in excellent agreement with the reference values obtained from independent simulations, for ΔT = ±15 K and ΔP = 100 bar (with respect to the simulated system). For larger systems, this range becomes smaller because the relative fluctuations of energy and volume become smaller. Without sufficient overlap, the performance of the method will become poor especially for nonlinear variations of the chemical potential.

We present a new molecular simulation code, Brick-CFCMC, for performing Monte Carlo simulations using state-of-the-art simulation techniques. The Continuous Fractional Component (CFC) method is implemented for simulations in the NVT/NPT ensembles, the Gibbs Ensemble, the Grand-Canonical Ensemble, and the Reaction Ensemble. Molecule transfers are facilitated by the use of fractional molecules which significantly improve the efficiency of the simulations. With the CFC method, one can obtain phase equilibria and properties such as chemical potentials and partial molar enthalpies/volumes directly from a single simulation. It is possible to combine trial moves from different ensembles. This enables simulations of phase equilibria in a system where also a chemical reaction takes place. We demonstrate the applicability of our software by investigating the esterification of methanol with acetic acid in a two-phase system.

In this paper, we review recent advances in the Continuous Fractional Component Monte Carlo (CFCMC) methodology and present a historic overview of the most important developments that have led to this method. The CFCMC method has gained attention for Monte Carlo simulations of adsorption at high loading, and phase and reaction equilibria of dense systems. It has recently been extended to reactive systems. The main features of the CFCMC method are: (1) Increased molecule exchange efficiency between different phases in single and multicomponent (reactive) systems, which improves the efficiency and accuracy of phase equilibria simulations at high densities; (2) Direct calculation of the chemical potential from a single simulation; (3) Direct calculation of partial molar properties from a single simulation. The developed simulation techniques are incorporated in the open-source molecular simulation software Brick-CFCMC.

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]

Deep eutectic solvents (DESs) are considered as green alternatives to room temperature ionic liquids (RTILs), due to their lower-cost synthesis and more environmentally friendly nature. In this work, Monte Carlo (MC) simulations have been used to compute the solubilities of CO₂, H₂S, CH₄, CO, H₂, and N₂ in choline chloride urea (ChClU) and choline chloride ethylene glycol (ChClEg) DESs. Due to the strong intermolecular interactions of DESs, leading to high viscosities, MC simulations present significant challenges with respect to system equilibration and solute molecule insertions. The Continuous Fractional Component Monte Carlo (CFCMC) method has been used with our open-source code, Brick-CFCMC, to improve molecule insertions and equilibration of the system, and directly compute the excess chemical potential and solubility (in terms of the Henry constant) of the gas molecules in the DESs. Pure DES properties, such as density and radial distribution functions (RDFs), were well reproduced by MC simulations. The solubilities of gases were, however, underestimated by the CFCMC simulations compared to available experimental data from literature. The order of solubilities of the different gases in ChClU at 328 K was obtained as H₂S > CO₂ > CH₄ > H₂ > CO > N₂, which reasonably agrees with experimental data from literature. The OPLS force field resulted in larger average Henry constants in ChClEg, compared to the GAFF force field, implying the better suitability of the GAFF force field for the calculations. Smaller ionic charge scaling factors were shown to increase the solubilities of the gases in the DESs, but result in lower densities. The differences between the computed Henry constants from MC simulations and experimental data from literature may be caused by the unsuitability of the used force field parameters of the DESs in combination with those of the solute gases. Nonetheless, experimental data from literature is scarce (except for CO₂) and in some cases contradictory, which makes the comparison with the computational results difficult.

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.

The CFCMC simulation methodology considers an expanded ensemble to solve the problem of low insertion/deletion acceptance probabilities in open ensembles. It allows for a direct calculation of the chemical potential by binning of the coupling parameter λ and using the probabilities p(λ = 0) and p(λ = 1), which require extrapolation. Here, we show that this extrapolation leads to systematic errors when the distribution p(λ) is steep. We propose an alternative binning scheme which improves the accuracy of computed chemical potentials. We also investigate the use of multiple fractional molecules needed in simulations of multiple components, and show that these fractional molecules are very weakly correlated and that calculations of chemical potentials are not affected. The statistics of Boltzmann averages in systems with multiple fractional molecules is shown to be poor. Good agreement is found between CFCMC averages (uncorrected for the bias) and Boltzmann averages when the number of fractional molecules is less than 1% of the total number of all molecules. We found that, in dense systems, biased averages have a smaller uncertainty compared to Boltzmann averages.

Hydrogen is one of the most popular alternatives for energy storage. Because of its low volumetric energy density, hydrogen should be compressed for practical storage and transportation purposes. Recently, electrochemical hydrogen compressors (EHCs) have been developed that are capable of compressing hydrogen up to P = 1000 bar, and have the potential of reducing compression costs to 3 kWh/kg. As EHC compressed hydrogen is saturated with water, the maximum water content in gaseous hydrogen should meet the fuel requirements issued by the International Organization for Standardization (ISO) when refuelling fuel cell electric vehicles. The ISO 14687-2:2012 standard has limited the water concentration in hydrogen gas to 5 μmol water per mol hydrogen fuel mixture. Knowledge on the vapor liquid equilibrium of H₂O-H₂ mixtures is crucial for designing a method to remove H₂O from compressed H₂. To the best of our knowledge, the only experimental high pressure data (P > 300 bar) for the H₂O-H₂ phase coexistence is from 1927 [J. Am. Chem. Soc., 1927, 49, 65-78]. In this paper, we have used molecular simulation and thermodynamic modeling to study the phase coexistence of the H₂O-H₂ system for temperatures between T = 283 K and T = 423 K and pressures between P = 10 bar and P = 1000 bar. It is shown that the Peng-Robinson equation of state and the Soave Redlich-Kwong equation of state with van der Waals mixing rules fail to accurately predict the equilibrium coexistence compositions of the liquid and gas phase, with or without fitted binary interaction parameters. We have shown that the solubility of water in compressed hydrogen is adequately predicted using force-field-based molecular simulations. The modeling of phase coexistence of H₂O-H₂ mixtures will be improved by using polarizable models for water. In the Supporting Information, we present a detailed overview of available experimental vapor-liquid equilibrium and solubility data for the H₂O-H₂ system at high pressures.

The combination of the TraPPE and OPLS/2016 force fields with five water models, TIP3P, SPC/E, OPC, TIP4P/2005 and TIP4P/EW was used to compute mixing enthalpies, excess chemical potentials and activity coefficients of water and methanol. Excess chemical potentials and activity coefficients were computed in an expanded version of the NPT ensemble. We found the best agreement between experimental data for all the computed properties of water–methanol mixtures for the TIP4P/2005-TraPPE force fields. The performance of the spherical cutoff methods in MC and MD simulations was compared to the Ewald summation. The radial distribution functions obtained from the Ewald summation and the Damped-Shifted Force (DSF) method were in excellent agreement. Numerical artifacts appeared at the cutoff radius when the original Wolf method was used to calculate the electrostatic interactions. The calculated excess mixing enthalpies, excess chemical potentials, and activity coefficients of water and methanol obtained from the Wolf method were in good agreement with the DSF method. Our simulation results show that the numerical artifacts of the original Wolf method have little effect for energy calculations in aqueous methanol mixtures.

An alternative method for calculating partial molar excess enthalpies and partial molar volumes of components in Monte Carlo (MC) simulations is developed. This method combines the original idea of Frenkel, Ciccotti, and co-workers with the recent continuous fractional component Monte Carlo (CFCMC) technique. The method is tested for a system of Lennard–Jones particles at different densities. As an example of a realistic system, partial molar properties of a [NH₃, N₂, H₂] mixture at chemical equilibrium are computed at different pressures ranging from P = 10 to 80 MPa. Results obtained from MC simulations are compared to those obtained from the PC-SAFT Equation of State (EoS) and the Peng–Robinson EoS. Excellent agreement is found between the results obtained from MC simulations and PC-SAFT EoS, and significant differences were found for PR EoS modelling. We find that the reaction is much more exothermic at higher pressures.

The applicability of the Wolf method for calculating electrostatic interactions is verified for simulating vapor-liquid equilibria of hydrogen sulfide, methanol, and carbon dioxide. Densities, chemical potentials, and critical properties are obtained with Monte Carlo simulations using the Continuous Fractional Component version of the Gibbs Ensemble. Saturated vapor pressures are obtained from NPT simulations. Excellent agreement is found between simulation results and data from literature (simulations using the Ewald summation). It is also shown how to choose the optimal parameters for the Wolf method. Even though the Wolf method requires a large simulation box in the gas phase, due to the lack of screening of electrostatics, one can consider the Wolf method as a suitable alternative to the Ewald summation in VLE calculations.