Multiple Free Energy Calculations from Single State Point Continuous Fractional Component Monte Carlo Simulation Using Umbrella Sampling.

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.

Introduction

Knowledge of chemical potentials is essential for understanding and studying phase behavior of multicomponent systems.[1,2] In Monte Carlo (MC) simulations, chemical potentials cannot be easily computed from a single simulation because chemical potentials cannot be expressed directly as a function of phase space variables (coordinates and momenta).[3−6] One of the most straightforward and commonly used methods to compute the chemical potentials is Widom’s test particle insertion (WTPI) method.[4,7,8] In this method, test particles are inserted at randomly selected positions inside the simulation box. The chemical potential is computed by sampling the interaction energies of the test particles with the other molecules in the system. In the NPT ensemble, the chemical potential of a certain component in a mixture is obtained fromin which ⟨···⟩ indicates an ensemble average, A is the component type, β = (kBT)−1, and kB is the Boltzmann constant. V and T are the volume and temperature of the system, respectively. ΔUA+ is the interaction energy of the test molecule of type A with the rest of the molecules. The partial molar volume and enthalpy of component type A can be written as derivatives of μA,[9] namely, h̅ = (∂(βμ)/∂β) and υ̅ = (∂μ/∂P).[9,10] In their pioneering work in the 80s of the past century, Frenkel et al. developed a method based on the WTPI method to compute partial molar properties.[11,12] These properties cannot be directly written as averages of phase space variables. Frenkel et al. showed that by random insertions of a test particle, partial molar properties (in the NPT ensemble) are obtained from[11]in which N denotes the number of molecules in the system and s denotes the scaled coordinates. In the case of an ideal gas, h̅Aex in eq becomes zero because there are no intermolecular interaction between the molecules. The derivation of eqs and 2b is shown in detail in our earlier work.[10] Partial molar properties are important for studying thermodynamic properties of multicomponent mixtures. Partial molar enthalpies also play a role in nonequilibrium thermodynamics, for example, in studying the heat flux in mixtures.[11] Recently, Josephson, Siepmann, and co-workers introduced a method to compute partial molar properties in the NPT version of the Gibbs ensemble using a method based on postprocessing fluctuations in energy, volume, and the number of molecules of each component.[13] It is also possible to compute partial molar properties from Kirkwood–Buff integrals.[14−20] In our recent work, we have outlined the importance of computing partial molar volumes υ̅ and partial molar enthalpies h̅ and provided an overview of different methods developed for computing partial molar properties in molecular simulation studies.[10] Similar to eq , efficient sampling of partial molar properties using eqs and 2b relies on the formation of cavities in the system, sufficiently large for a test particle. The probability of formation of such cavities is directly related to the density of the system. The limitations of insertion/deletion type methods are well-known for dense systems and/or systems with complex molecules.[4,10,21−26] At high density, test particle insertions result in overlaps with existing molecules in the system, leading to a very large potential energy ΔUA+, that is, a very small Boltzmann factor as shown in eqs , 2a and 2b. Because of increased computational power since the early 1950s,[27−29] more complex systems have been investigated in molecular simulation studies of phase equilibria.[30−33] The inherent limitations of the WTPI method in free energy calculations in complex systems have led to the development of other simulation methods.[7,8,34]

In the previous century, alternative methods were introduced to circumvent the sampling problems of the WTPI method, such as the overlapping distribution method (ODM) by Bennett[35] or combining particle insertion and deletion methods by Shing and Gubbins.[36−38] The ODM is based on the free energy difference between two systems. This computation is not efficient if the systems significantly differ in configuration space.[7,25] Other Widom-type methods such as the Theodorou deletion method[22] and Widom’s test particle deletion method[21] suffer from the same sampling difficulties associated with the WTPI method for complex and dense systems.[39] Lyubartsev et al. developed expanded ensembles for accurate free energy calculations of complex systems.[3] This method takes advantage of gradual insertion/deletion of molecules which allows the system to gradually adjust to the new configuration in which a new molecule is added.[40−46] Inspired by this, Shi and Maginn developed the continuous fractional component MC (CFCMC) method to facilitate molecule exchanges in open ensemble simulation of dense phases.[47] In CFCMC simulations, a fractional molecule is added to the ensemble. During the simulation, the interactions of this molecule are adjusted with a coupling parameter (λ). The fractional molecule acts as ideal gas if λ = 0 and interacts fully with other molecules if λ = 1. The application of the CFCMC method for systems ranging from Lennard-Jones (LJ) particles to larger molecules such as hydrocarbons is well-established, see refs (48) and (49). Although the method by Shi and Maginn facilitates molecule exchanges in dense phases, postprocessing is required to obtain the chemical potentials. Recently, Vlugt and co-workers combined the CFCMC method with free energy calculations in ref (2). In this method (hereafter referred to as the CFCMC method), the chemical potential can be directly computed by sampling the probability distribution p(λ). Consider a multicomponent mixture of mono-atomic components, expanded with a fractional molecule. In the NPT ensemble, the ensemble partition function can be written as[10,50]in which UfracA is the interaction energy of the fractional molecule. Depending on the value of λ, UfracA can vary between “zero” and “fully scaled”.[1,1,47,51] The chemical potential of component A in CFCNPT equals[50]in which ρA is the number density of A. To make the logarithm in eq dimensionless, an arbitrary reference density ρ0 is selected. The term μ0 contains intramolecular contributions to the chemical potential. p(λA = 1) and p(λA = 0) are the probabilities of λA when λA = 1 and λA = 0, respectively. The chemical potential of eq is split into an ideal gas part and excess part. It is shown in ref (10) that the computed chemical potentials obtained using eqs and 4 are the same by definition. In the CFCNPT ensemble, the partial molar excess enthalpy and partial molar volume are equal[10]in which ⟨H⟩ and ⟨V⟩ are ensemble averages of the total enthalpy and total volume computed at λA = 1 and ⟨H/V⟩ and ⟨1/V⟩ are ensemble averages computed at λA = 0. In ref (10), we showed that eqs , 2b, 5a, and 5b are identical and yield identical results. Because the intramolecular contribution is not included in QCFC of eq , the excess partial molar enthalpy computed using eq is always taken with respect to the ideal gas. Because the WTPI method suffers at high densities, the application of eqs , 2a, and 2b is rather limited, in sharp contrast to eqs , 5a, and 5b. In our recent publications, the advantages of using the CFCMC method for free energy calculations are shown for various systems.[1,2,10,52−55] To further use the CFCMC method, we introduce a formulation, based on the idea of umbrella sampling,[56] allowing one to compute the chemical potentials (eq ) for multiple state points from a single state point simulation.

Umbrella sampling is a well-known method developed by Torrie and Valleau[56] from which the free energy difference between the system of interest and a reference system can be obtained.[7,34,56] By introducing a biasing function , the ensemble average ⟨A⟩ in the canonical ensemble is calculated from an ensemble where the configuration space is sampled proportional to in which ⟨···⟩Π denotes an ensemble average in the ensemble Π. The biasing function is only a function of the configuration space. The derivation is provided in ref (34). In this paper, we show that by combining eq with CFCMC, the chemical potentials can be accurately estimated for an appreciable temperature and pressure range from a single state simulation in the CFCNPT ensemble. To estimate the chemical potential at a temperature or pressure different from that of the simulation, sufficient overlap is required between the configuration space of the two systems.[7,34] This depends also on the system size. The relative density and energy fluctuations become smaller with the increase in system size, which may reduce the overlap between the configuration spaces of two systems. Therefore, we investigated the limitations of combining umbrella sampling with the CFCMC method when estimating chemical potentials for various temperature and pressure intervals. The combination of umbrella sampling with the CFCMC method offers the following advantages: (1) Accurate estimates of the chemical potential as a function of pressure or temperature are computed from a single simulation. (2) Partial molar properties are obtained directly from a single simulation. By definition, partial molar properties can be obtained by numerically evaluating the expressions h̅ = (∂(βμ)/∂β) and υ̅ = (∂μ/∂P), using the estimated chemical potentials at different pressures and temperatures.[9,10] (3) Partial molar enthalpies obtained from umbrella sampling can be used as an independent check for the CFCMC method in eqs and 5b, both for programming bugs or an independent check whether the phase space is sufficiently sampled.

In Section , the combination of umbrella sampling and free energy calculations in CFCMC simulations is explained. In Section , we provide a detailed overview of the systems considered in our simulations. In Section , our simulation results are presented. In this section, the chemical potentials for a binary LJ color mixture are estimated at multiple temperatures and pressures by performing single state simulations in the CFCNPT ensemble. The results from umbrella sampling are in excellent agreement with the ones obtained from independent simulations. In addition, partial molar properties are estimated by numerically evaluating h̅ = (∂(βμ)/∂β) and υ̅ = (∂μ/∂P), based on the results from umbrella sampling. The computed partial molar properties for the LJ mixture are in excellent agreement with the results from the CFCMC method (eqs and 5b) and the WTPI method (eqs and 2b). As an example of a molecular system, we applied our method to different mixtures of water and methanol at standard conditions (T = 298 K and P = 1 bar) and compared the results with those obtained from independent CFCNPT simulations. Accurate estimates of the chemical potential of water and methanol are obtained using umbrella sampling for a temperature difference ΔT = ±10 K, for N = 470 molecules. Excellent agreement is observed between partial molar properties of water and methanol obtained from umbrella sampling and those obtained using eqs and 5b. We also applied our method to a mixture of ammonia, nitrogen, and hydrogen at T = 573 K and P = 800 bar at chemical equilibrium and compared the partial molar properties obtained from our method to those obtained from the CFCMC method. The limitations of the method are tested by decreasing overlap in configuration space for pure methanol. In Section , our conclusions are presented.

Theory

Similar to eq , we combine umbrella sampling with the CFCMC method to estimate the probability distribution p(λ) at (T, P★) while performing a simulation at (T★, P★) in the CFCNPT ensemble. This distribution is calculated usingin which H = U + PV, U is the internal energy of the system, and c is a normalization constant. In the ensemble, at β*, the value of the scaling parameter equals λ′. In the Supporting Information, one can find the derivation of eq . To estimate the excess part of the chemical potential in eq at β, the values p(λ = 1)|β and p(λ = 0)|β are obtained using eq . Similar to eq , it can be shown that the number density can be calculated at T* usingin which ρ|β is the number density at β. The number of fractional molecules is not included in computing number densities.[1] If the number of fractional molecules is much smaller compared to that of the whole molecules (less than 1%), Boltzmann averages are not affected.[50] By combining eqs , 7, and 8, the chemical potential at (T, P★) is estimated while performing a single simulation at (T★, P★). In a similar manner, the distribution p(λ) at (T★, P) can be estimated while performing a simulation at (T★, P★), keeping the number of molecules constant. This distribution is obtained by introducing a bias to the ensemble average

The derivation of eq can be found in the Supporting Information. By computing the values p(λ = 0)| and p(λ = 1)| based on eq , one can obtain the excess chemical potential at pressure P while performing a simulation at pressure P★. For the ideal part of the chemical potential in eq , the number density is estimated at P using umbrella sampling

By combining eqs , 9, and 10, the chemical potential at (T★, P) is obtained while performing a single state simulation at (T★, P★). The derivation of eq is very similar to the derivation of eq .

Simulation Details

Besides thermalization trial moves in the CFCNPT ensemble, three additional sorts of trial moves are defined to facilitate gradual insertions/deletions of the fractional molecule: (1) changes in λ: randomly increasing or decreasing λ while the positions and orientations of all molecules are fixed. For the fractional molecule, the LJ interactions are scaled using[57,58]

It is possible to choose different alchemical pathways for uLJ. A detailed study of different alchemical pathways can be found in refs (57−59). For details on scaling, how the electrostatic interactions are scaled, the reader is referred to refs (1, 52, 55). (2) Identity changes: an attempt is made to transform the fractional molecule into a whole molecule and simultaneously transform a randomly selected whole molecule into a fractional molecule. (3) Random re-insertions: the fractional molecule is re-inserted at a randomly selected position. Except for the fractional molecule, orientations and positions of other molecules do not change. The acceptance rules for these trial moves are derived based on metropolis importance sampling.[1,27] The trial move (2) has a high acceptance probability at high values of λ, and the trial move (3) has a high acceptance probability at low values of λ. Therefore, it is efficient to combine trial moves (2) and (3) and use a single hybrid trial move.[1,10] It is shown in ref (10) that combining the last two trial moves leads to higher acceptance ratios of trial moves (2) and (3). To flatten the λ space, prior to performing production runs, a weight function W(λ) was computed using the Wang–Landau (WL) algorithm.[60,61]

To demonstrate the feasibility of our method, MC simulations were carried out to compute the chemical potentials of an LJ color mixture (50–50%) in the CFCNPT ensemble, containing N = 200 molecules. All simulations were carried out at T* = 2 and pressures between P* = 0.1 and P* = 8 (the symbol * denotes the reduced units). For the LJ system, eqs and 8 were used to compute the chemical potential at temperatures between T* = 1.82 and T* = 2.22 and a pressure of P* = 6 from a single simulation at (T* = 2, P* = 6). Similarly, eqs and 10 were used to compute the chemical potentials at T* = 2 and pressures between P* = 5.95 and P* = 6.05 in the same simulation. To compute the reference values for the chemical potentials, independent simulations were performed for each temperature and pressure. The estimated chemical potentials obtained from umbrella sampling were used to numerically evaluate the partial molar properties using the expressions h̅ = (∂(βμ)/∂β) and υ̅ = (∂μ/∂P). As a reference, the partial molar properties were also computed using the WTPI method using eqs and 2b and the CFCMC method using eqs and 5b. All LJ interactions were truncated and shifted at 2.5σ, and the thermal wavelength Λ was set to unity. For every simulation of the LJ system in the CFCNPT ensemble, 105 cycles were used for equilibration and 107 cycles were used for production. The number of trial moves in each cycle is equal to the total number of molecules in the system. The following probabilities were set for selecting trial moves: 33% translations, 33% λ changes, 33% hybrid trial moves, and 1% volume changes.

To test this method for a more complex system, aqueous methanol mixtures with mole fractions of methanol between xMeOH = 0.2 and xMeOH = 0.8 were simulated in the CFCNPT ensemble. All MC simulations were performed using our in-house code. We verified that the same results were obtained for various systems[1,52] when the RASPA software package[62,63] was used. For water, the TIP4P/2005[64] and for methanol, the TraPPE[65] force fields were used. Sites with a partial charge but no LJ parameters are protected by adding an LJ site with ϵ/kB = 1 K and σ = 1 Å. A comparison between different combinations of force fields for predicting thermodynamic properties of aqueous methanol mixtures is provided in our earlier work.[55] In all simulations of water–methanol mixtures, a fractional molecule of methanol and a fractional molecule of water were present. All molecules were modeled as rigid objects, and pairwise interaction potentials consisted only of LJ and Coulomb potentials. A cutoff radius of 12 Å for both LJ and Coulomb potentials was used. Lorentz–Berthelot was used for cross interactions between methanol and water molecules, and analytic tail corrections were applied.[7,8] The damped shifted force method[66] is used to handle the electrostatic interactions of water and methanol systems. For details on using these cutoff-based methods, the reader is referred to refs (52, 55, 66−69). The interactions of the fractional molecule are scaled in the following way: starting from λ = 1, the electrostatic interactions are scaled to zero and thereafter the LJ interactions are scaled to zero. This circumvents configurations with charge overlaps between the fractional molecule and other molecules leading to numerical instabilities in the simulation.[57−59,70−72] Chemical potentials of water and methanol between T = 288 K and T = 308 K were estimated by performing umbrella sampling from a single simulation of the water–methanol mixture at T = 298 K and P = 1 bar. As a reference, the chemical potentials of water and methanol were computed independently at every pressure and temperature. Subsequently, partial molar properties of water and methanol were obtained from the estimated chemical potentials using umbrella sampling. The results were compared to the partial molar properties obtained using eqs and 5b.[10] For every simulation of the water–methanol system in the CFCNPT ensemble, 105 cycles were used for equilibration and 6 × 106 cycles for production. The following probabilities were set for selecting trial moves: 35% translations, 30% rotations, 17% λ changes, 17% hybrid trial moves, and 1% volume changes. The WL algorithm[60,61] was used independently for water and methanol to flatten the probability of λ.[60,61] In ref (50), we showed that the correlation between scaling parameters of different fractional molecules is weak and independent of the weight function. It is therefore computationally advantageous to separate the multidimensional weight function into a series of one dimensional weight functions: W(λ1, λ2) ≈ W(λ1) + W(λ2). This is due to the fact that filling one-dimensional histograms is more efficient than filling two separate multidimensional histograms.[50] To investigate the limitations of umbrella sampling in CFCMC simulations, umbrella sampling was applied to estimate the chemical potentials of pure methanol for a wide pressure and temperature range. The results were compared to the results obtained from independent simulations of pure methanol in the CFCNPT ensemble, at every temperature and pressure.

To compute partial molar properties of ammonia, nitrogen, and hydrogen in their mixture at T = 573 K and P = 800 bar, simulations were performed in the CFCNPT ensemble.[10] Similar to our previous work, separate simulations were performed in which only a single fractional molecule was present. This was repeated for the three components, leading to three independent simulations. Force field and simulation details are provided in ref (10). The starting mixture composition (NNH = 407, NN = 7, and NH = 20) was obtained from simulations of the Haber–Bosch process in the reaction ensemble as described in ref (1).

Results and Discussion

A simulation in the CFCNPT ensemble was performed for an LJ system at T* = 2 and P* = 6. Using umbrella sampling, the chemical potentials were computed for temperatures ranging from T* = 1.82 to T* = 2.22 and pressures ranging from P* = 5.95 to P* = 6.05. The values of the chemical potential were compared to those obtained from independent CFCNPT simulations for all temperatures and pressures. The results are shown in Table . Excellent agreement is observed between the computed chemical potentials obtained using umbrella sampling from a single simulation at (T* = 2, P* = 6) and the chemical potentials calculated from independent simulations at P* = 6 and temperatures between T* = 1.82 and T* = 2.22. In this work, the error bars (uncertainties) are obtained by computing the standard deviation of the mean from five independent simulations. It should be noted that the uncertainty of the computed chemical potential associated with umbrella sampling increases with the increase in the temperature difference |ΔT*|. Because the overlap between the configuration spaces decreases with the increase in temperature or pressure difference, it impairs the sampling. The overlap between the configuration spaces also decreases with the increase in system size, as the relative fluctuations of energy and volume become smaller for larger systems. In Table b, the estimated chemical potentials from umbrella sampling from a single simulation at (T* = 2, P* = 6) are compared to the chemical potentials obtained from independent simulations at T* = 2 and pressures between P* = 5.95 and P* = 6.05. Excellent agreement is observed between the computed chemical potentials estimated using umbrella sampling and the results obtained from independent simulations for the entire pressure range. From the results presented in Table , it can be observed that accurate estimates of the chemical potentials are obtained from a single simulation combining the CFCMC method with umbrella sampling. Although the statistical uncertainties of the estimated chemical potentials obtained from umbrella sampling are larger compared to those obtained from independent simulations, the differences in uncertainties are not significant.

Table 1

Chemical Potentials for a Binary (50–50%) LJ Mixture Obtained from Umbrella Sampling (Eqs and 9) and Independent Simulations in the CFCNPT Ensemble (Eq )a

To obtain estimates of the chemical potentials at (a) different temperatures and (b) pressures, umbrella sampling is performed at T* = 2 and P* = 6 in a single simulation. Boltzmann averages at T* = 2 and P* = 6 are highlighted in gray. Numbers in brackets indicate uncertainties.

The chemical potentials obtained from umbrella sampling are used to compute the partial molar properties for the LJ mixture, by numerically differentiating (∂(βμ)/∂β) and (∂μ/∂P). For instance, the chemical potentials shown in Table are used to plot βμ–β and μ–P for T* = 2 and P* = 6, as shown in Figure . The data points connected with a line as shown in Figure indicate that the data points were obtained from a single simulation combining CFCMC and umbrella sampling, and the individual data points were obtained from independent CFCNPT simulations. Calculating the slopes as shown in Figure , at T* = 2 and P* = 6, leads to the values of the partial molar volumes and enthalpies. A central difference scheme with high order approximation was used to compute the partial molar enthalpies as shown in Figure a and partial molar volumes as shown in Figure b. The partial molar excess enthalpy was computed by subtracting the ideal gas contribution from the total partial molar enthalpy. It is shown in ref (10) that the ideal part of the partial molar enthalpy of an LJ particle, h̅id equals 5/(2β). It is noteworthy that the contribution of the thermal wavelength Λ in h̅id equals 3/(2β).[10] Because Λ is set to unity in our simulations, the contribution of the thermal wavelength cancels out when numerically evaluating (∂(βμ)/∂β). Therefore, the partial molar excess enthalpy is obtained from umbrella sampling of the chemical potential using h̅ex = (∂(βμ)/∂β) – 1/β. The partial molar properties for the LJ system obtained from umbrella sampling and the CFCMC method (eqs and 5b) and the WTPI method (eqs and 2b) are provided in Table . Excellent agreement is observed between the partial molar properties from CFCNPT simulations and the results obtained from umbrella sampling. Similar uncertainties are observed between the two methods. For the LJ system, umbrella sampling works equally well compared to the CFCMC method. In addition, accurate estimates of the chemical potential at different pressures and temperatures are obtained without any extra computational power or post processing.

Figure 1

Plots showing (a) β*μ–β* and (b) μ–P* for an LJ binary color mixture consisting of 200 molecules. Downward-pointing triangles indicate the results obtained from independent CFCNPT simulations at pressures between P* = 5.94 and P* = 6.05 and temperatures between T* = 1.82 and T* = 2.22. Circles indicate the results obtained by performing umbrella sampling from a single simulation in the CFCNPT ensemble at T* = 2 and P* = 6. Lines indicate that the data are obtained from a single simulation. Error bars are smaller than the symbol sizes. Raw data are listed in Table .

Table 2

Densities, Chemical Potentials, Partial Molar Excess Enthalpies, and Volumes at T* = 2 and Pressures between P* = 0.1 and P* = 8 Computed for a Binary LJ Color Mixture (200 Molecules)a

			CFCNPT		umbrella sampling		WTPI method
P*	⟨ρ*⟩	μA	h̅Aex	υ̅A	h̅Aex	υ̅A	h̅Aex	υ̅A
0.1	0.052	–7.460(7)	–0.44(4)	18.6(3)	–0.44(5)	19(2)	–0.361(5)	19.29(4)
2	0.584	–1.075(8)	–1.6(1)	1.72(3)	–1.6(1)	1.72(3)	–1.64(4)	1.70(1)
4	0.722	1.957(5)	–0.06(5)	1.39(1)	–0.06(5)	1.39(1)	–0.3(2)	1.36(5)
6	0.800	4.581(9)	1.7(2)	1.26(2)	1.7(2)	1.26(2)	1.6(8)	1.23(5)
8	0.856	7.001(6)	3.2(1)	1.17(1)	3.2(1)	1.17(1)	4(1)	1.2(1)

Three different methods are used at each P*, in order to compare the results for the partial molar properties: the WTPI method (eqs and 2b) in the NPT ensemble,[11,12] the CFCMC method (in the CFCNPT ensemble, eqs and 5b),[10] and umbrella sampling (eqs and 9). The chemical potentials are calculated using eq from independent simulations. Numbers in brackets indicate uncertainties.

Other realistic systems with complex intermolecular interactions considered here are different aqueous methanol mixtures with mole fractions of methanol between xMeOH = 0.2 and xMeOH = 1. These systems were simulated in the CFCNPT ensemble. For a mixture composition xMeOH = 0.8, umbrella sampling was used to estimate the chemical potentials of water and methanol at temperatures between T = 288 K and T = 308 K while running a single CFCNPT simulation at P = 1 bar and T = 298 K. In addition, independent simulations in the CFCNPT simulations were performed to obtain the chemical potentials as a reference. The results are shown in Table . Accurate estimates of the chemical potentials of water and methanol are obtained between T = 288 K and T = 308 K, from a single CFCNPT simulation T = 298 K (ΔT = ±10 K). The relative differences between the estimates of the chemical potentials and the chemical potentials obtained from independent simulations are well below 1%. For each temperature, the differences between the absolute values of the chemical potentials as shown in Table are significantly smaller than 1 kcal/mol = 4.184 kJ/mol, which is typically considered as a benchmark in the computational chemistry literature.[73] The results show that umbrella sampling can provide an accurate estimate of the chemical potentials of water and methanol for an appreciable temperature range (around ΔT = ±10 relative to the simulation temperature). The raw data shown in Table were used to plot βμ as a function of β for the water–methanol mixture (N = 470), xMeOH = 0.8, at temperatures between T = 288 K and T = 308 K. The results are shown in Figure . The lines indicate data obtained from a single simulation. It can be seen in Figure that excellent agreement is observed between the results obtained from umbrella sampling and the reference values from independent simulations in CFCNPT at every temperature. The error bars associated with the estimated chemical potentials for water and methanol increase with increase in |ΔT| relative to T = 298 K. This indicates that umbrella sampling becomes inefficient and inaccurate for a larger temperature difference. It should be emphasized that this temperature range is not a priori known and is system- and system size-dependent. We explain later in this paper how a reasonable temperature or pressure range can be selected for umbrella sampling of the chemical potentials.

Table 3

Chemical Potentials [kJ/mol] of the Water–Methanol Mixtures (xMeOH = 0.8) Obtained from Umbrella Sampling and Independent Simulations in the CFCNPT Ensemble at Temperatures between T = 288 K and T = 308 Ka

The values of the chemical potentials are estimated using umbrella sampling at T = 298 K and P = 1 bar. The corresponding Boltzmann averages are highlighted in gray. Numbers in brackets indicate uncertainties.

Figure 2

Plot showing βμ–β for a water–methanol mixture, xMeOH = 0.8. To compute the chemical potentials, independent simulations are performed at temperatures between T = 288 K and T = 308 K. Downward-pointing triangles and upward-pointing triangles denote data for water and methanol from independent CFCNPT simulations. Asterisks connected by a line indicate data obtained using umbrella sampling from a single simulation at T = 298 K and P = 1 bar. Error bars are smaller than the symbol sizes. Raw data are listed in Table .

The βμ–μ and μ–P plots for the water–methanol mixture at xMeOH = 0.8 are shown in Figure . The lines indicate that the data points are obtained from a single simulation. To compute the partial molar properties of water and methanol at T = 298 K, a central difference scheme with high order approximation is used to numerically evaluate (∂βμ/∂β) and (∂μ/∂P). For each component type, the partial molar excess enthalpy of water and methanol is then computed by subtracting the ideal gas part. For different mixture compositions, partial molar enthalpies of methanol and water are computed using umbrella sampling and used as a reference. The results are compared to those obtained from the CFCMC method. The results are shown in Table for the mixture with mole fractions of methanol between xMeOH = 0.2 and xMeOH = 1. Because the density of methanol mixtures at T = 298 K is high, the WTPI method is not used here. Excellent agreement is found between both methods. It is observed both for the LJ system and water–methanol mixture that the accuracy and uncertainty of the partial molar properties from umbrella sampling and the CFCMC method are similar. It is also observed that the uncertainties associated with partial molar properties are an order of magnitude larger compared to the uncertainties of the chemical potentials. This is due to the potential drawback of eqs and 5b when subtracting two large numbers with a (relatively) small difference.[10] This may induce larger error bars compared to the chemical potential. A similar potential drawback is observed when using umbrella sampling to compute partial molar properties. Numerically computing the derivatives (∂(βμ)/∂β) and (∂μ/∂P) also involves subtracting two numbers (or several numbers) over a relatively small temperature or pressure difference. Therefore, accurate estimates of the chemical potential are needed to obtain accurate values for the partial molar properties. Based on the results, it is clear that the overall performance of both methods is very similar in terms of accuracy and precision.

Table 4

Chemical Potentials [kJ/mol], Partial Molar Excess Enthalpies [kJ/mol], and Partial Molar Volumes [cm3/mol] for Water–Methanol Mixtures at Different Compositions at T = 298 K and P = 1 bar, Using Simulations in the CFCNPT Ensemble (Eqs and 5b) and Umbrella Sampling (Eqs and 9)a

				CFCNPT		umbrella sampling
i	xi	ρi	μitot	h̅iex	υ̅i	h̅iex	υ̅i
MeOH	0.2	289.4(3)	–34.0(1)	–40(10)	40(10)	–40(10)	40(10)
H2O	0.8	650.9(3)	–38.9(1)	–47(7)	17(14)	–47(7)	17(14)
MeOH	0.4	485.7(6)	–32.7(2)	–42(9)	39(9)	–42(8)	39(9)
H2O	0.6	409.6(5)	–39.3(1)	–47(7)	16(5)	–47(7)	16(5)
MeOH	0.6	620.7(4)	–32.2(1)	–41(5)	40(5)	–41(5)	40(5)
water	0.4	233.0(2)	–39.7(1)	–49(9)	11(7)	–49(9)	11(7)
MeOH	0.8	713.7(6)	–31.5(1)	–41(6)	45(5)	–41(6)	45(5)
H2O	0.2	100.3(1)	–41.4(3)	–47(12)	21(16)	–47(12)	21(16)
MeOH	1	777.0(9)	–31.1(1)	–42(3)	37(5)	–42(3)	37(5)

Umbrella sampling is performed at T = 298 K and P = 1 bar for all mixture compositions. x denotes the mole fraction of the component i. The chemical potentials are calculated from eq using independent simulations. Numbers in brackets indicate uncertainties.

To investigate the limitations of umbrella sampling in free energy calculations, a wide temperature and pressure range is selected for estimating the chemical potential of pure methanol from a single CFCNPT simulation. The chemical potentials of pure methanol in the temperature range of T = 266 K and T = 340 K and the pressure range of P = 1 bar to P = 1001 bar were computed using eqs –8 from a single simulation at T = 298 K and P = 1 bar. It should be noted that a wide temperature and pressure range is only selected to test the limitations of the method. The temperature and pressure range should be selected to ensure sufficient overlap in configuration space and energy between the systems. As a reference, independent CFCNPT simulations of methanol were performed at each temperature and pressure to compute the corresponding chemical potentials. The results are shown in Table . It can be seen that the estimated chemical potentials for ΔT = ±15 K, relative to T = 298 K, are in excellent agreement with the chemical potentials obtained from independent simulations. It is clear that the sampling becomes more difficult when the temperature difference increases. The uncertainties of the estimated chemical potentials are an order of magnitude higher for |ΔT| > 15 K. This sampling difficulty can be explained based on the overlap between the energy or configuration space of a system at two different temperatures or pressures. In Figure a, the probability distribution of enthalpy per molecule p(h) for pure methanol at different temperatures is shown (h = H/N, N = 410). For the distribution p(h) at T = 266 K, no overlap is observed with the distribution p(h) at T = 298 K, and the method fails to estimate the chemical potential at T = 266 K, as can be seen in Table a. At T = 283 K, sufficient overlap is observed with the distribution p(h) at T = 298 K, and the estimated chemical potential is in excellent agreement with the reference value. It should be noted that the uncertainties of the chemical potentials from independent simulations are always smaller. However, for ΔT = ±15 K, the uncertainties of the estimated chemical potentials is small. In Figure a, it is also shown that the overlaps between distributions p(h) for T = 298 K and T = 320 K become very small. This leads to large uncertainties in the estimated values for the chemical potentials. It is shown in Figure b that for P < 100 bar, the estimated chemical potentials of methanol are in excellent agreement with the chemical potentials obtained from independent simulations. For pressures ranging between 100 and 500 bar, the uncertainties of the results from umbrella sampling are up to 2 orders of magnitude larger than those obtained from independent simulations. The estimated values of the chemical potentials are however within chemical accuracy (1 kcal/mol).[73] In Figure b, the probability distributions of volume per molecule p(υ) for pure methanol at different pressures are shown (υ = V/N, N = 410). Sufficient overlap is observed in Figure b for the distributions p(υ) at P = 1 bar and P = 100 bar. This is expected because the compressibility of liquid methanol is very low at room temperature. Therefore, it is possible to estimate the chemical potentials of methanol accurately for any pressure between P = 1 bar and P = 100 bar, see Table b and Figure b. For P > 500 bar, the uncertainty of the estimated chemical potentials increases significantly (3 orders of magnitude larger compared to independent simulations) and the method starts to break down. This can be explained by examining the overlap between the distributions p(υ) at high pressures and P = 1 bar. In Figure b, no significant overlap is observed between the distribution p(υ) at P = 1 and the distributions at P > 500. Under these conditions, the excess chemical potentials of methanol computed from umbrella sampling deviate significantly from the excess chemical potentials computed from independent simulations. Here, we illustrate the sampling issue of the excess chemical potential at high pressures by plotting p(λ) as a function of pressure, computed from umbrella sampling at T = 298 K and P = 1 bar. For pure methanol (N = 410), it is clearly observed in Figure that the computed p(λ) shows scatter for pressures significantly different from that of the simulation. This clearly illustrates that the method breaks down when the pressure difference becomes large, and the statistics for the excess chemical potential become very poor. We observed that for liquid methanol (N = 410), the computed 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. The pressure and temperature range for accurate estimation of the chemical potentials or other thermodynamic properties may differ from one system to another. However, plotting distributions p(h) and p(υ) for a wide pressure and temperature range can readily visualize at what range umbrella sampling can be applied to obtain accurate results. As shown in Figure , an investigation of the overlap between energies, volumes at different temperatures, and pressures can easily indicate the boundaries at which the method starts to fail. For systems that lack sufficient overlap (i.e., Δβ or ΔP is too large), it is expected that the performance of the method is poor. If one is interested in nonlinear variations of μ (e.g., higher order derivatives of μ), our method may not work well because of a lack of overlap. To compute partial molar properties, selecting a narrow region (a small temperature or pressure range) still allows for a numerical evaluation of h̅ = (∂(βμ)/∂β) and υ̅ = (∂μ/∂P). Therefore, partial molar properties can be computed from umbrella sampling without selecting a wide temperature or pressure range. As an example of a strong hydrogen bond-forming system, we considered the Haber–Bosch process (N2 + 3 H2 ⇌ 2NH3). Umbrella sampling was used to compute partial molar properties of the mixture at T = 573 K and P = 800 bar, at chemical equilibrium. From our earlier work, the equilibrium composition of the mixture is known.[1,10] Independent simulations are performed at T = 573 K and P = 800 bar in the CFCNPT ensemble. In each simulation, the fractional molecule of only one component is present. Partial molar properties obtained from umbrella sampling (eqs and 9) and the CFCMC method (eqs and 5b) are compared in Table . Excellent agreement is observed between the results obtained from umbrella sampling and the CFCMC method.

Table 5

Chemical Potentials [kJ/mol] of Pure Methanol Obtained from Umbrella Sampling (Eqs –8) and Independent Simulations in the CFCNPT Ensemble (Eq )a

Umbrella sampling is performed in the CFCNPT ensemble of methanol at T = 298 K and P = 1 bar. Boltzmann averages obtained from umbrella sampling are highlighted in gray. Numbers in brackets indicate uncertainties.

Figure 3

(a) Probability distribution of the enthalpy per molecule of methanol; h = H/N at P = 1 bar and different temperatures: T = 266 K (green), T = 285 K (teal), T = 298 K (black), T = 320 K (orange), and T = 340 K (red). (b) Probability distribution of the volume per molecule of methanol; υ = V/N at T = 298 K and different pressures: P = 1 bar (black), P = 100 bar (teal), P = 500 bar (orange), and P = 1000 bar (red). The number of molecules for pure methanol in the liquid phase is N = 410.

Figure 4

Probability distribution of λ for pure methanol obtained from a CFCNPT simulation at T = 298 K and P = 1 bar. Equation is used to compute p(λ) for pressures up to P = 600 bar from a single simulation. The red lines indicate the Boltzmann distribution p(λ) at P = 1 bar and the distribution p(λ) computed for P = 600 bar.

Table 6

Partial Molar Enthalpies [kJ/mol] and Partial Molar Volumes [cm3/mol] of Ammonia, Nitrogen, and Hydrogen at T = 573 K and P = 800 bara

	CFCNPT		umbrella sampling
	h̅	υ̅	h̅	υ̅
NH3	–10(1)	50(7)	–10(1)	50(7)
N2	11(1)	108(6)	11(1)	108(6)
H2	12(1)	94(6)	12(1)	94(6)

The composition of the mixture is obtained from simulations of ammonia synthesis in the reaction ensemble, as explained in our earlier work.[1] Partial molar properties are computed in the CFCNPT ensemble using eqs and 5b and umbrella sampling using eqs and 9. Numbers in brackets indicate uncertainties.

Conclusions

We introduced an alternative method to obtain accurate estimates of the excess chemical potential of a component for a wide temperature and pressure range from a single simulation. This method combines umbrella sampling and the CFCMC technique. Using the values of the estimated chemical potentials, the partial molar enthalpies and volumes of a component are obtained by numerically evaluating the derivatives (∂(βμ)/∂β) and (∂μ/∂P), respectively. This method does not have the disadvantages of the WTPI method. As a proof of concept, the values of the chemical potential for a binary LJ mixture were estimated using umbrella sampling from a single simulation in the CFCNPT ensemble at T* = 2 and P* = 6. For a temperature range between T* = 1.82 and T* = 2.22 and pressure range between P* = 5.95 and P* = 6.05, excellent agreement was observed between the estimated chemical potentials and those obtained from independent CFCNPT simulations of the LJ system, at each temperature and pressure. Partial molar properties obtained from umbrella sampling were in excellent agreement with the CFCMC method introduced in ref (10) and the original method of Frenkel et al.[11,12] We observed that the accuracy and precision of the averages obtained from the CFCMC method in ref (10) and umbrella sampling are very similar. To apply our method to a system with complex intermolecular interactions, we considered aqueous mixtures of methanol with different compositions. We applied our method to estimate the chemical potentials of methanol and water for a temperature range between T = 288 K and T = 308 K from a single simulation at T = 298 K and P = 1 bar. Excellent agreement was found between the results obtained from umbrella sampling and those obtained from independent simulations in the CFCNPT ensemble, with relative differences well below 1%. As an example of a strong hydrogen bond-forming system, our method was applied to a mixture of ammonia, nitrogen, and hydrogen at chemical equilibrium at T = 573 K and P = 800 bar. It was observed that partial molar properties of ammonia, nitrogen, and hydrogen obtained from umbrella sampling and the CFCMC method are in excellent agreement. We investigated the limitations of our method for liquid methanol (N = 410 molecules). It was observed that for the temperature difference ΔT = ±15 K, very accurate estimates of chemical potential at different temperatures were obtained from umbrella sampling from a single CFCMC simulation. In addition, it was found that for a pressure difference ΔP = 100 bar, accurate estimates of chemical potential at different pressures were obtained from umbrella sampling from a single simulation. Lack of sufficient overlap (ΔP or Δβ) between different states may result in a poor performance of the method. Based on the results, it can be concluded that combining umbrella sampling with CFCMC provides a powerful tool for accurate estimation of chemical potentials for an appreciable temperature and pressure range and computation of partial molar properties. This can potentially reduce the number of simulations of a system for free energy calculations for a specific temperature and/or pressure range.