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.

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.

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.

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.

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

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.