Molecular Simulation of Vapor-Liquid Equilibria Using the Wolf Method for Electrostatic Interactions.

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.

Introduction

Vapor–liquid equilibria (VLE) are of great interest in the chemical industry, in particular for separation processes.[1,2] As an alternative to experiments, molecular simulation is a useful tool for obtaining VLE data, especially for multicomponent systems.[3−5] Molecular simulation critically relies on a force field and the geometry of the molecules as input. The force field describes the interactions between atoms and molecules. The past decades, there have been many simulation studies on VLE using force fields that employ a combination of a Lennard-Jones potential (LJ) and an electrostatic potential.[6−10] Different methods and force fields for describing interactions exist, such as density functional theory[11,12] and polarizable force fields.[13−15] The LJ potential is rather short-ranged, and calculating the energy is therefore straightforward, that is, by truncating the LJ potential and using analytic tail corrections.[16] However, the electrostatic potential is long-ranged and one has to take extra care when calculating this energy contribution because of periodic boundary conditions and the interaction of the molecule with its mirror images. The Ewald summation[17] is a commonly used method for calculating electrostatics. It calculates the electrostatic energy by dividing the potential in a short-range potential that can be calculated directly and a long-range potential which requires a Fourier transform of the charge density. The Ewald summation is accurate for a wide variety of systems and the electrostatic energy is well-defined.[17,18] A disadvantage of the Ewald summation and its variants is that it is computationally expensive because a Fourier transform is involved.[19,20] Other methods for dealing with electrostatics exist such as the reaction field method,[21,22] the particle–particle and particle–mesh algorithm[23] (a variant to the Ewald summation), and the Wolf method.[24] Here we focus on the latter. The past few years, the Wolf method has become quite popular for simulating dense fluids due its simplicity and efficiency.[25] The Wolf method makes use of the (strong) screening of electrostatic interactions in dense systems. All interactions are pairwise and there is no Fourier transform of the charge density involved. Calculations involve two parameters: the damping parameter α and the cutoff radius Rc. Unlike the Ewald method, there is no clear criterion for choosing these parameters. The effectivity of the Wolf method has already been demonstrated in many applications.[26−30] Those simulations take place in dense liquids where there is a strong screening of the electrostatic interactions. From earlier studies, we also know that the Wolf method is not very efficient in zeolites because of less effective screening due to voids in the zeolite topology.[31] As a consequence, a larger cutoff radius for electrostatic interactions is needed in this case, making the Wolf method less efficient than one would expect. It is also unclear if and how the Wolf method influences vapor–liquid equilibria. Here, we will investigate the accuracy of the Wolf method in simulating VLE for simple compounds: hydrogen sulfide, methanol, and carbon dioxide. These compounds were chosen because of their relevance in industry and applications and because accurate force fields are available to describe these compounds. Monte Carlo simulations are performed in the Continuous Fractional Component (CFC) version of the Gibbs Ensemble (GE). In this ensemble, the system is expanded by a so-called fractional molecule that has scaled interactions with the surrounding particles.[32] The advantages of this ensemble are that molecule insertions are much more efficient than in the conventional GE, and one can directly obtain the chemical potential of the system, even at high densities. Simulating in the GE makes it possible to choose the parameters for the Wolf method (α and Rc) in both boxes (vapor and liquid) independently of each other. This is important because the Wolf method is based on the screening of charges which is dependent on the density of the system.

This article is organized as follows. First, the methodology is described and the (scaled) potentials are defined as well as the expression for obtaining chemical potentials. Second, the simulation details, such as the input and force fields are discussed. Also the procedure on how to choose the optimal Wolf parameters is explained. Third, the simulation results for density, chemical potential, vapor pressure, and critical properties are summarized. We finish with our concluding remarks about the efficiency and accuracy of the Wolf method in VLE simulations.

Methodology

We simulate vapor–liquid coexistence by Monte Carlo simulations in the GE using the CFCMC method by Poursaeidesfahani and co-workers.[32] Although the multiple histogram reweighting method[33,34] is a very efficient method for obtaining VLE, MC simulations in the GE are suitable for obtaining VLE with a relatively small amount of molecules not too close to the critical point.[35] In the GE two boxes are present: one for the liquid phase and one for the vapor phase. A simulation involves displacements and rotations of molecules, volume changes of the boxes, and molecule transfers between the boxes. In the CFC version of the GE, one extra molecule per molecule type is added to the system which we call the fractional molecule. Interactions of this molecule are scaled by a parameter λ which ranges from 0 (no interactions with surrounding molecules) to 1 (full interactions with surrounding molecules). The fractional molecule can be in either of the simulation boxes. Besides the standard trial moves for thermalization: translations, rotations, and volume changes, there are three additional trial moves to facilitate molecule transfer between the boxes: (1) a change in the value of the fractional parameter λ; (2) a swap move which transfers the fractional molecule from one box to the other; (3) an identity change move where the fractional molecule turns into a whole one, and a molecule in the other box turns into a fractional molecule. For the acceptance rules regarding the trial moves that facilitate molecule transfers we refer to ref (32). The CFC version of the GE is more efficient than the standard GE because of higher acceptance probabilities for the molecule transfer. It has been shown that the presence of a fractional molecule has a negligible effect on the thermodynamic properties of the system.[36] The CFCMC simulations require a (two-dimensional) weightfunction W(λ, i) (i = 1,2 indicates the box where the fractional molecule is) to flatten the λ-probability distribution so that all values of λ are equally likely to occur and the system does not get stuck at a certain value (or range) of λ. Furthermore, W is used to make the fractional molecule be equally likely to be found in one of the simulation boxes. Molecules in the simulations are considered as rigid and interactions are described by a combination of electrostatic and LJ potentials.

Electrostatics

Electrostatic interactions are calculated using the Wolf method:where q is the partial charge of atom i, erfc is the complementary error function, α is a damping parameter, r = |r⃗ – r⃗| is the distance between two atoms i and j, Rc is the cutoff radius, and N is the total number of atoms in the system. We would like to point out that a modification of the Wolf method exists[25] which has a continuous transition in the electrostatic force and can therefore be used in Molecular Dynamics simulations. Since we are here only interested in the thermodynamics and not the dynamics of the system, we use the Wolf method in its original form. The values of α and Rc can be chosen independent in the two boxes in the GE. As shown later, for calculating the correct electrostatic energy, the value of α needs to be larger in the box with the liquid phase than in the box with the gas phase, and the value of Rc is smaller in the liquid phase than the gas phase. This is related to the fact that in the liquid phase the effective screening is larger and the interaction is more short-ranged. The electrostatic interaction of the fractional molecule is scaled by substituting r → r + A(1 – λ)2 in eq and substituting q → λq in eq and eq for the atoms of the fractional molecule. The term A(1 – λ)2 prevents singularities for very small values of r. The value of A is chosen as Å.

Lennard-Jones Interaction

For the LJ interactions (truncated at R):where ϵ and σ are the LJ-parameters between atom i and j, analytic tail corrections are used:[16]where the sum ranges over all atom types in the system, N is the number of atoms of type i (excluding the fractional molecules), and V is the volume of the simulation box. The factor 1/2 accounts for double counting interactions. The LJ interactions of the fractional molecules in the CFCMC method are scaled according to[37,38]and the fractional molecule contributes to the energy tail correction by substituting N → N + λ in eq . This ensures that for λ = 0, the fractional molecule does not contribute to the tail correction, and for λ = 1 the contribution is equal to that of a normal molecule in the system.

Chemical Potential

In the CFCMC method the total chemical potential of a pure-component system in box i can be calculated according to[32]where V, N, and p(λ) are the volume, the number of (whole) particles, and the probability distribution of λ in box i, respectively. We use the notation p(λ ↑ 1) and p(λ ↓ 0) for the limit p(λ) approaching λ = 1 from below and p(λ) approaching λ = 0 from above. As p(λ) can be quite steep, especially near λ = 1, a quadratic extrapolation is used to obtain those limits. For convenience, the thermal wavelength Λ is set to 1 Å for all systems.

Simulation Details

The VLE of three different pure compounds are simulated: hydrogen sulfide, methanol, and carbon dioxide. The force field parameters[3,39,40] used to describe the interaction between molecules are listed in Table and Lorentz–Berthelot rules are used.[16] The molecules are treated as rigid and bond lengths and angles are listed in Table . Periodic boundary conditions are used.

Table 1

Force Field Parameters of Hydrogen Sulfide,[39] Methanol,[40] and Carbon Dioxide.[3] M is a Dummy Site, CH3 Is Described as United Atom

	atom	ϵ/kB/K	σ/Å	q/e
H2S	S	122.0	3.60	0
	H	50.0	2.50	0.21
	M	0	0	–0.42
CH3OH	CH3	110.45	3.6499	0.1546
	O	97.775	3.1659	–0.6544
	H	0	0	0.4998
CO2	C	27.0	2.800	0.700
	O	79.0	3.050	–0.350

Table 2

Geometries of Hydrogen Sulfide,[39] Methanol,[40] and Carbon Dioxide[3]

		atoms	bond length/angle
H2S	bond	S–H	1.34 Å
	bond	S–M	0.30 Å
	angle	H–S–H	92°
	angle	H–S-M	46°
CH3OH	bond	CH3–O	1.43 Å
	bond	O–H	0.945 Å
	angle	CH3–O–H	108.5°
CO2	bond	C–O	1.16 Å
	angle	O–C–O	180°

A MC cycle in the simulations consist of Ntot MC steps, where Ntot is the total number of molecules in the whole system. In each MC step, a trial move is selected at random with the following fixed probabilities: 35% translations, 30% rotations, 1% volume changes, 17% λ changes, 8.5% swaps of the fractional molecule, and 8.5% identity changes of the fractional molecule. Ensemble averages are updated after every MC cycle. For the fractional molecule, starting from λ = 0, first the LJ interaction is increased gradually from no interaction to full interaction with surrounding molecules and after that the electrostatic interaction is increased from no interaction to full interaction at λ = 1.[41]

Each simulation starts with 5 × 103 MC-cycles for equilibrating, where only translation and rotation moves are performed. Next, there are 5 × 104 cycles for initializing and further equilibrating of the system, now using all trial moves. In this phase, the weightfunction W(λ, i) is being constructed using the Wang–Landau algorithm.[42] During the initializing and equilibrating phases, the maximum displacement, rotation, and volume change are modified to achieve an acceptance ratio of 50% for those trial moves. Finally, there are 5 × 104 production cycles where ensemble averages are taken and the λ-probability distribution is sampled for which 100 bins are used for storage.

Multiple simulations for different compounds and different Wolf parameters at different temperatures are performed. For all simulations a cutoff radius of 14 Å for LJ interaction is used in both boxes with tail corrections. The Wolf parameters (α and Rc) are chosen for each simulation as follows. First, a short NVT simulation is run at a density close to equilibrium (estimated from literature[39,43,44]) above the critical temperature (so that no phase separation occurs in the simulation box). Second, for the final configuration of this NVT simulation the electrostatic energy is calculated with the Ewald summation as well as for many different Wolf parameters. A plot comparing the electrostatic energy calculated with the Ewald summation (which we consider as the exact solution) and the Wolf method for different parameters is made. Figures and 2 show typical differences in the electrostatic energy for the liquid and vapor phase of methanol, respectively. Figure also clearly shows the effect of a damping parameter that is too small: for α → 0 and large Rc the lack of screening in the cutoff spheres result in large energy differences. From those plots, the optimal values for α and Rc can be determined; that is, the parameters that give an accurate result compared to the Ewald summation, choosing Rc as small as possible. We will refer to these values as the optimal parameter set. It is possible to take multiple configurations and use the averages of the energy differences to obtain the parameter sets. However, it is sufficient to take only one configuration as can be seen for example in Figure where, for α large enough, the energy differences for different Rc values converge. We also verified that the optimal parameter set is the same as when we would have taken a configuration from a GE simulation after equilibrating and follow the same procedure. As choosing the best values for α and Rc can be quite some work we also run simulations for different (easier) choices of the parameters to show the dependence of the results on the values of the Wolf parameters. We consider the following sets of Wolf parameters:The second set is chosen because it should give accurate results in the liquid phase of each compound (according to the optimal set). The third set is chosen with a slightly smaller α such that it should be more accurate for the gas phase and is a value in the range that is typically chosen.[25,26,31] The fourth set is chosen such that it should be even more accurate for the gas phase. In the optimal parameter set for very low densities an α of 0 is found with a very large cutoff radius. This corresponds to calculating the electrostatic energy directly, ignoring any interaction between a molecule and its mirror images.

Figure 1

Relative difference in electrostatic energy between the Wolf method and Ewald summation for different values of Rc as a function of α. The parameters for the Ewald summation are chosen such that a relative precision of 10–6 is achieved. The energy is calculated for methanol at a (typical liquid) density of 692 kg/m3 at 600 K. The optimal value of α is in the range from 0.10 Å–1 to 0.14 Å–1.

Figure 2

Relative difference in electrostatic energy between the Wolf method and Ewald summation for different values of Rc as a function of α. The parameters for the Ewald summation are chosen such that a relative precision of 10–6 is achieved. The energy is calculated for methanol at a (typical gas) density of 2.66 kg/m3 at 600 K. The optimal value of α is in the range from 0 to 0.03 Å–1.

optimal (as described in the text)

α = 0.12 Å–1 and Rc = 14 Å (16 Å for CO2) in both boxes

α = 0.10 Å–1 in both boxes and the same Rc as in the optimal set

α = 0.06 Å–1 in both boxes and the same Rc as in the optimal set

Each simulation starts with different initial configurations (how the molecules are distributed over the two boxes and the size of the boxes) which can be found in the Supporting Information. In the SI also the optimal Wolf parameter sets can be found for each system.

After obtaining the densities at equilibrium, critical properties and the saturated vapor pressures Pvap are determined. The critical properties are calculated from the VLE-curves using the method described by Dinpajooh and co-workers[35] as well as in the book by Frenkel and Smit.[45] The saturated vapor pressures are determined in the following way. We estimate the pressure (for example from the ideal gas law) and set up multiple NPT simulations at different pressures P in a range around the estimated one. From each simulation we calculate the density and from those results a (P,ρ)-diagram can be constructed. Pvap can then be determined (by interpolating) at the equilibrium vapor density ρv.

All simulations were performed with an in-house code, optimized for VLE simulations. It was verified to produce the same results as the RASPA software package[46,47] with the Ewald summation. The reported uncertainties in the results are standard deviations obtained from results from six independent simulations.

Results

The densities at different temperatures can be found in Figure for hydrogen sulfide, Figure for methanol, and Figure for carbon dioxide. Our simulation results are in excellent agreement with the comparing data (simulations using the Ewald method) for almost all Wolf parameter sets. Only at high densities does the density deviate a bit for the set where α = 0.06 Å–1 in both boxes. This is caused by the fact that at high density there is a large effective screening of charges. This means that either a larger value of α should be used (see also Figure ). Closer to the critical temperatures small deviations in the results are found.

Figure 3

Density–temperature plot for the vapor–liquid equilibrium of hydrogen sulfide. Different symbols indicate different parameters for the Wolf method, compared to simulation results from Shah[39] using the Ewald summation (solid line). Tabulated data together with the uncertainties can be found in the Supporting Information (Table S4).

Figure 4

Density–temperature plot for the vapor–liquid equilibrium of methanol. Different symbols indicate different parameters for the Wolf method, compared to experimental results from Goodwin[43] (solid line) and simulations from Gonzalez–Salgado[40] using the Ewald summation (dashed line). Tabulated data together with the uncertainties can be found in the Supporting Information (Table S5).

Figure 5

Density–temperature plot for the vapor–liquid equilibrium of carbon dioxide. Different symbols indicate different parameters for the Wolf method, compared to experimental results and equation of state from Duschek[50] (solid line) and simulations from the NIST database[44] using the Ewald summation (dashed line). Tabulated data together with the uncertainties can be found in the Supporting Information (Table S6).

The total chemical potential of hydrogen sulfide, methanol, and carbon dioxide can be found in Figures , 7, and 8, respectively, as calculated from eq . From this data it can immediately be seen that the vapor and liquid are in chemical equilibrium because the chemical potential is equal in both boxes. For methanol there is a slight difference between the chemical potential in the liquid and the gas phase, but this is still within the error bars. Furthermore, for methanol at T = 180 K and T = 230 K, no chemical potential could be calculated directly because of the very low density of the gas phase. Although it would be possible to calculate the chemical potential for the gas phase directly by considering it as an ideal gas, we do not do that here because of the large error bars for the density.

Figure 6

Total chemical potential as a function of temperature of hydrogen sulfide. Values are obtained with eq where Λ = 1 Å. Tabulated data together with the uncertainties can be found in the Supporting Information (Table S7).

Figure 7

Total chemical potential as a function of temperature of methanol. Values are obtained with eq where Λ = 1 Å. Tabulated data together with the uncertainties can be found in the Supporting Information (Table S8).

Figure 8

Total chemical potential as a function of temperature of carbon dioxide. Values are obtained with eq where Λ = 1 Å. Tabulated data together with the uncertainties can be found in the Supporting Information (Table S9).

Critical temperatures and densities are determined from the VLE-curves. The obtained critical points of hydrogen sulfide, methanol, and carbon dioxide are listed in Table . For hydrogen sulfide and carbon dioxide acceptable critical points are obtained. Only for the Wolf parameter set with α = 0.06 Å–1 in both boxes a larger deviation is found. This can again be explained by the fact that this parameter set fails to give an accurate description of the VLE (Figure ). For methanol, the critical temperatures are higher than expected and the error bars are relatively high compared to that for hydrogen sulfide and carbon dioxide. This is most likely caused by the difference in the shape of the VLE-curve, compared to that of hydrogen sulfide and carbon dioxide, and the method used to extract the critical point.

Table 3

Critical Temperatures Tc and Critical Densities ρc for Different Compounds and Different Wolf Parameter Setsa

	H2S		CH3OH		CO2
	Tc	ρc	Tc	ρc	Tc	ρc
parameter set	K	kg·m–3	K	kg·m–3	K	kg·m–3
1	377(4)	345(8)	536(16)	257(14)	307(2)	463(10)
2	374(5)	346(10)	533(19)	261(13)	307(3)	463(7)
3	375(2)	347(5)	537(21)	260(19)	308(2)	464(10)
4	378(8)	347(17)	546(10)	262(6)	320(3)	457(5)
ref (sim)	374.5(6)	361(2)	521.5	272	306.2	464.9
ref (exp)	373.4(2)	347(4)	513(1)	273(3)	304.18(2)	468(1)

Critical properties are determined from the VLE curves with the method described by Dinpajooh and co-workers.[35] Parameter set 1 is the optimal parameter set. Parameter set 2 has α = 0.12 Å–1 and Rc = 14 Å (16 Å for CO2) in both boxes. Parameter set 3 has α = 0.10 Å–1 in both boxes and the same Rc as in the optimal set. Parameter set 4 has α = 0.06 Å–1 in both boxes and the same Rc as in the optimal set. The last rows show reference (simulation and experimental) data for hydrogen sulfide,[39,48] methanol,[40,44] and carbon dioxide.[3,49] Critical properties are determined using data in the following temperature ranges: hydrogen sulfide 260–320 K, methanol 360–480 K, and carbon dioxide 220–280 K. The numbers between brackets denote the uncertainty in the last digit.

Finally, from NPT simulations the vapor pressures at equilibrium are determined. Only the optimal Wolf parameter sets were used. A Clausius–Clapeyron plot (Figure ) summarizes the resulting vapor pressures. We were not able to determine the vapor pressure from a NPT simulation at T = 200 K because no density corresponding to ρv was obtained. Also, at T = 360 K we were not able to obtain the pressure because the temperature is close to the critical point. For the same reason we were not able to obtain the vapor pressure for carbon dioxide at T = 300 K.

Figure 9

Clausius–Clapeyron plots of the saturated vapor pressure as a function of the inverse temperature. The solid lines show the results obtained from NPT simulations using the Wolf method for electrostatic interactions with the optimal parameters, different colors indicate the different compounds. Dashed lines show simulation data from the NIST database,[44] Shah,[39] and Gonzalez-Salgado.[40] Tabulated data together with the uncertainties can be found in the Supporting Information (Table S10).

The optimal set of Wolf parameters produces results in agreement with VLE data from literature and can therefore be considered as an accurate alternative to the Ewald method. Also note that the other parameter sets produce acceptable results as well, especially Set 2. This is most likely caused by the fact that it is more important to describe interactions in the liquid phase more accurately than it is in the vapor phase (which is close to behaving like an ideal gas). This is a useful result if one does not have data available for determining the optimal parameters α and Rc but still wants to perform simulations using the Wolf method. Simulations using the Wolf method are found to be at least twice as fast in terms of CPU time compared to the Ewald method.

Conclusion

We tested the applicability of the Wolf method for electrostatic interactions in MC simulations. By performing simulations in the GE with the CFCMC method, we were able to obtain accurate densities, chemical potentials, critical points, and vapor pressures at VLE for different compounds and different Wolf parameter sets. We showed that the damping parameter α in the liquid phase should be chosen larger than α in the vapor phase, corresponding to a more effective screening of charges in the liquid phase. The cutoff radius Rc can be chosen smaller in the liquid phase than Rc in the gas phase, because of the higher value of α in the liquid. Moreover, we showed that a simple estimation of the Wolf parameters can already produce accurate results for VLE.