Brick-CFCMC: Open Source Software for Monte Carlo Simulations of Phase and Reaction Equilibria Using the Continuous Fractional Component Method.

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.

Introduction

Knowledge of phase and reaction equilibria is essential for the design and operation of chemical processes.[1−4] In the past decades, molecular simulation has become an important tool for predicting these equilibria based on the interactions between molecules: the so-called force field.[5−8] While there are many software packages available for molecular dynamics simulations (e.g., Gromacs,[9,10] LAMMPS,[11,12] NAMD[13]), the number of open-source Monte Carlo codes for studying phase and reaction equilibria is limited[14] (e.g., Towhee,[15] Cassandra,[16] RASPA,[14,17] FEASST,[18] GOMC[19]), and rarely gradual insertions and deletions of molecules are considered. Here, we present Brick-CFCMC (hereafter referred to as Brick), a new open-source software package for force field-based Monte Carlo simulations in various ensembles (e.g., NVT/NPT, the Gibbs Ensemble,[20,21] the Reaction Ensemble,[22−24] and the Grand-Canonical/Osmotic Ensemble[25−28]).

Having sufficient molecule exchanges is often the most critical part of Monte Carlo simulations in open ensembles. Important features of Brick are as follows: (1) Molecule exchanges are facilitated by the use of fractional molecules,[29] which significantly improves insertion and deletion of molecules and allows for direct calculation of chemical potentials[30] (and their derivatives[31]). (2) Flexibility of molecules is taken into account by bond-bending, torsion, and intramolecular nonbonded interactions. (3) Intermolecular and intramolecular interactions are described by a combination of Lennard-Jones and electrostatic interactions. (4) Both the Ewald method[32] and Wolf method[33] can be used for electrostatic interactions and (5) smart Monte Carlo trial moves,[34,35] enabling collective displacements and rotations. Single- and multicomponent systems can be used. Brick has been used in some of our recent simulation work, e.g., the study of the ammonia synthesis reaction,[36] the computation of partial molar properties,[31] the computation of solubility of water in high-pressure hydrogen,[37] the study of vapor–liquid equilibria of xylene mixtures,[38,39] and the computation of thermodynamic properties of water–methanol mixtures.[40] Brick is open source and can be downloaded from https://gitlab.com/ETh_TU_Delft/Brick-CFCMC.

Software

Brick is written in Fortran 90/95 and works most conveniently in workflow that involves a terminal environment. It can be compiled using the GNU Fortran Compiler[41] or the Intel Fortran Compiler.[42] The package includes the main source code for performing molecular simulations as well as tools that help to create input for the simulations. For more information, we refer the reader to the Supporting Information and the Brick manual.[43]

The Continuous Fractional Component (CFC) method[29−31,36,44] is implemented in this software. This method significantly improves the efficiency of particle transfers by introducing so-called fractional molecules to the system. To each of these fractional molecules, a continuous fractional parameter, λ ∈ [0, 1], is assigned that regulates the strength of interactions of the fractional molecule with the surrounding molecules. When λ = 0, no interactions are present which means that this molecule is treated as an ideal gas molecule. This enables efficient molecule insertions and deletions. At λ = 1, the fractional molecule has the same interactions as the other molecules not being a fractional molecule, so-called whole molecules. This enables trial moves that facilitate the change of identity of the fractional molecules with whole molecules, which in turn facilitate molecule exchanges. Introducing the parameter λ adds an additional degree of freedom to the system. It has been shown that adding fractional molecules does not affect the results of the simulation.[45,46] From the probability distribution of λ, the CFC method allows for direct calculation of chemical potentials (μ),[30,31,36,47] fugacity coefficients (φ),[31,36] partial molar enthalpies (h̅),[31,47] and partial molar volumes (v̅).[31,47] To prevent the system from getting stuck at certain values of λ, biasing of λ is necessary. This is achieved by using a weight function (W(λ)) that can be obtained via the Wang–Landau algorithm[48] or an iterative scheme.

The following ensembles are available for performing simulations: (1) the NVT and NPT ensembles, with and without the CFC method, (2) the Gibbs Ensemble,[21,30] for simulations of phase equilibria, at constant total volume or at constant pressure, (3) the Reaction Ensemble,[22,24,36,44] for simulations of reaction equilibria, at constant volume or at constant pressure, and (4) the Grand-Canonical Ensemble,[29] at constant volume or at constant pressure. The latter is also referred to as the Osmotic Ensemble and is only meaningful for multicomponent systems in which one extensive variable is fixed.[49] In Brick, it is possible to combine ensembles to simulate a combined phase and reaction equilibrium.

For electrostatic interactions, it is possible to use the Ewald method, the Wolf method,[33] or the damped-shifted version of the Wolf method.[50] For liquids in which electrostatic interactions are well screened, the Wolf method is computationally more efficient than the Ewald method while leading to the same accuracy.[50,51]

Other features of Brick are the calculation of the partition function of isolated molecules from spectroscopic data[52] or from QM simulations (e.g., Gaussian[53]) required for simulations in the Reaction Ensemble, calculation of fugacities in multicomponent systems using the Peng–Robinson Equation of State,[54] calculation of radial distribution functions, and Widom’s test-particle insertion method.[31,55]

Case Study: Esterification of Acetic Acid with Methanol

To demonstrate features of Brick, we simulate the esterification of acetic acid with methanol in the liquid phase[56,57]Two phases form in this reaction due to the polar nature of acetic acid, methanol and water, and the nonpolar methyl acetate. Chemical reaction equilibria can be simulated in the Reaction Ensemble, where molecules are converted from reactants into products (and vice versa) using Monte Carlo trial moves such that the reaction equilibrium is obtained.[22,24,36,44] Phase equilibria can be simulated in the Gibbs Ensemble, where molecules are being transferred between two simulation boxes such that the phase equilibrium is obtained.[21,30] The Gibbs Ensemble is combined with the Reaction Ensemble for simulating the esterification. This is illustrated in Figure . Details are provided in the Supporting Information.

Figure 1

Schematic representation of the combination of the Gibbs Ensemble with the Reaction Ensemble for the esterification of methanol with acetic acid, CH3OH + CH3COOH ⇄ CH3COOCH3 + H2O. A total of 900 molecules is distributed over the two simulation boxes. Four fractional molecules, one for each component, are added to the system to facilitate molecule transfers between the two simulation boxes. These fractional molecules can be in either simulation box and can change from one to the other simulation box by the Gibbs Ensemble swap and identity change trial moves.[30] Fractional molecules of either reactants or products are added to each simulation box to facilitate the chemical reaction. These fractional molecules remain in the same simulation box during the simulation and can be converted from reactants into products (or vice versa) by the Reaction Ensemble swap and identity change trial moves.[36] This means that in total, at all times during the simulation, we have eight fractional molecules (four for achieving phase equilibrium and two in each simulation box to achieve reaction equilibrium). This is less than 1% of the total number of molecules so that the fractional molecules do not affect the simulation results.[46] Figure created with iRASPA.[58]

The equilibrium composition obtained from the simulations are listed in Table . We observe a clear phase separation. One phase contains mostly water, while the other phase contains mostly methyl acetate. We refer to these two phases as the water-rich phase and ester-rich phase, respectively. The chemical potential per mole can be written as the sum of three parts[31,36]where R is the gas constant; T is the temperature; q is the ideal gas partition function, excluding the translational part, of an isolated molecule of component i; V0 is an arbitrary reference volume (here set to 1 Å3); Λ is the thermal wavelength; ⟨ρ⟩ is the average the number density; ρ0 is an arbitrary reference number density (here set to 1 Å–3), and p(λ = 1) and p(λ = 0) are the probabilities that the fractional parameter λ takes the values 1 and 0, respectively. The first term on the right-hand side only depends on the temperature and is denoted by μ°. The second term depends on the number density (and is together with the first term referred to as the ideal gas contribution to the chemical potential). The third term is the excess chemical potential. Chemical potentials that were obtained from simulations in the NPT ensemble are listed in Table . We observe that for each component the chemical potentials in both phases are equal. This shows that the system is at phase equilibrium. To verify that the system achieved reaction equilibrium, we add the chemical potentials (eq ) of the reactants (CH3OH + CH3COOH) and of the products (CH3COOCH3 + H2O) in both phases (Supporting Information). We conclude that we achieved reaction equilibrium in each simulation box, phase equilibrium between the simulation boxes, and that flexibility of the molecules does not affect the results.

Table 1

Compositions and Chemical Potentials at Reaction and Phase Equilibrium of the Esterification of Methanol with Acetic Acid at T = 343 K and P = 1 bara

	Component	xiI	xiII	μiI – μi° (kJ·mol–1)	μiII – μi° (kJ·mol–1)	ai
Rigid	CH3OH	0.135(9)	0.072(6)	–36.0(8)	–35.8(5)	0.15(2)
	CH3COOH	0.06(2)	0.151(7)	–38.1(9)	–38.7(7)	0.21(5)
	CH3COOCH3	0.14(2)	0.67(2)	–29.6(10)	–30.3(5)	0.9(2)
	H2O	0.67(2)	0.11(2)	–37.0(4)	–36.6(7)	1.4(2)
Flexible	CH3OH	0.150(7)	0.073(8)	–35.7(8)	–36.2(6)	0.14(3)
	CH3COOH	0.06(2)	0.17(2)	–39.0(10)	–39.1(6)	0.17(5)
	CH3COOCH3	0.120(9)	0.69(3)	–30.1(8)	–30.1(8)	0.9(2)
	H2O	0.66(2)	0.08(3)	–37.1(4)	–37.9(5)	1.1(2)

Results for a system where all molecules are rigid and a system where molecules are flexible (i.e. bond bending and torsion are taken into account). The superscripts I and II indicate the water-rich and ester-rich phases, respectively, x is the mole fraction of component i, μ is the chemical potential (eq ), μ° = −RT ln(qV0/Λ3) is the contribution to the chemical potential due to the internal degrees of freedom (listed in the Supporting Information), and a is the thermodynamic activity. The number between brackets indicates the uncertainty (one standard deviation) in the last digit.

The equilibrium constant of the reaction can be computed fromWe obtain ln K = 3.8 ± 0.4 for the system with rigid molecules, and ln K = 3.7 ± 0.5 for the system with flexibile molecules. The uncertainty is defined as the standard deviation. The computed equilibrium constants are in good agreement with experimental data[59] where values for ln K are found to be in the range from 3.25 to 3.41. The equilibrium constant would be ln KIG = 2.48 if the system was treated as an ideal gas, showing that the medium has a large effect on the reaction equilibrium.

Conclusion

We present a new open-source software package, called Brick-CFCMC, for molecular simulations of phase and reaction equilibria using state-of-the-art Monte Carlo simulation techniques. Chemical potentials, fugacity coefficients, and partial molar properties are directly calculated from single simulations. We point out the applicability of Brick-CFCMC to industrial relevant processes by a study of the esterification of methanol with acetic acid. For this system, we calculated the equilibrium composition, chemical potentials, thermodynamic activities, and activity coefficients. We conclude that there is no significant difference in treating the molecules as rigid or flexible for this system.