LETTER pubs.acs.org/JPCL
Using an Analytic Equation of State to Obtain Quantitative Solubilities of CO2 by Molecular Simulation Christian S. Schacht,† Thijs J. H. Vlugt,† and Joachim Gross*,‡ † ‡
Process & Energy Laboratory, Delft University of Technology, Leeghwaterstraat 44, 2628 CA Delft, The Netherlands Institute of Thermodynamics and Thermal Process Engineering, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany
bS Supporting Information ABSTRACT: The design of solvents for CO2 capture could significantly profit from the use of appropriate simulation tools to accurately determine solute-solvent properties. Force-field-based molecular simulations can be applied to obtain accurate predictions of thermophysical properties of pure components. However, force fields often fail to correctly describe interactions in binary and multicomponent systems, especially for mixtures of polar and nonpolar components. Here, we present a method to quickly generate optimized force field parameters by using a statistically based equation of state. The corrected parameters are used in computer simulations to calculate the solubility of carbon dioxide in liquid n-alkanes over a wide range of temperatures. The predicted values of the Henry coefficient are in excellent agreement with experimental data for all systems investigated. The presented formalism can be applied to any force field and property of interest if a statistically based model is available to describe this property. SECTION: Statistical Mechanics, Thermodynamics, Medium Effects
T
he European Union is committed to a 20% reduction of their carbon dioxide emissions (as compared to the 1999 level) by the year 2020.1 Emission control and, thus, carbon dioxide capture and storage (CCS) are considered key elements for reaching these objectives.2,3 The development of new solvents for carbon capture and knowledge about their affinity for carbon dioxide is crucial to establish cost-efficient CCS systems. Force-field-based molecular simulation is a powerful tool for the evaluation and prediction of solvent properties4 as its successful application significantly reduces the number of time-consuming experiments. In addition, these simulations enable us to obtain a molecular understanding of the key factors that lead to a desired solubility. In these simulations, we assume a model (often called force field) to describe the interactions between the atoms in the system. Many force fields are designed to be transferable, that is, the parameters for a given interaction site can supposedly be transferred to various molecules. Results of molecular simulation critically rely on the quality of the used force field. Unfortunately, direct application of commonly used force fields often result in inaccurate predictions of solubilities in complex systems4-6 as these force fields are usually parametrized only from pure component vapor-liquid equilibria. This seriously hinders the successful application of molecular simulation to design solvents with certain properties. In this study, we show that the statistically based PCP-SAFT equation of state7 can be used to overcome this limitation by efficiently generating transferable force field parameters for r 2011 American Chemical Society
mixtures. We derive a direct relation between the PCP-SAFT parameters and the force field parameters. We then use the PCPSAFT equation of state to model the outcome of molecular simulations. Consequently, the optimization of the force field can be performed within the PCP-SAFT framework. The resulting parameters can be applied in molecular simulations and lead to excellent agreement with experimental data. It is important to note that evaluating PCP-SAFT is many orders of magnitude faster than performing a molecular simulation. This leads to a major reduction of computational efforts when compared to conventional methods for the optimization of force fields5,8 and enables an efficient application of computer simulations for the molecular design of solvents. In this study, we apply this approach to predict the solubility of carbon dioxide in liquid n-alkanes, for which the well-known TraPPE force field9,10 results in significant deviations from experiments. By calibrating the dispersive energy for the alkane-CO2 interactions using the PCP-SAFT equation of state, we obtain quantitative predictions of carbon dioxide solubilities in liquid n-alkanes over a wide temperature range. The PCP-SAFT equation of state provides an analytical expression for the Helmholtz energy of a system of interacting molecules. Molecules are described as chains of tangent spherical segments of the same size that interact with a 12-6 Lennard-Jones Received: December 14, 2010 Accepted: January 19, 2011 Published: February 02, 2011 393
dx.doi.org/10.1021/jz101688h | J. Phys. Chem. Lett. 2011, 2, 393–396
The Journal of Physical Chemistry Letters
LETTER
(united) atoms of type R in molecule type i. The factor ψij is introduced to eq 3 to correct for model inaccuracies and approximations that entered the derivation of this equation (see the Supporting Information). The factor ψij is usually close to unity. Once the constant factor ψij is determined, we can use PCPSAFT to model the outcome of a force-field-based molecular simulation for the observable Ωsim(p) as a function of the force field parameters p that we wish to optimize. The optimization of the force field parameters p can consequently be performed using PCP-SAFT only. For binary systems of carbon dioxide and n-alkanes, the terms describing the cross-component interaction can be identified from eq 3 as
(LJ) potential. Despite its simplicity, the essential characteristics of repulsive and dispersive interactions of nonspherical molecules are captured by this model. The molecular model can be refined further by adding associating interactions11 of short-range or multipolar potentials.12,13 Within PCP-SAFT, a nonassociating, nonpolar molecule is defined by the number of segments m̂i, the segment size parameter σ̂i, and a dispersive energy parameter ε̂i. In sharp contrast to PCP-SAFT, the force-field-based approach used in molecular simulations evaluates the total potential energy of the system as a function of the positions of all atoms in the system. Molecular simulations are used to compute thermodynamic properties of such a system. The well-known TraPPE united atom force field9 describes CHx groups as single interaction sites. The nonbonded interactions between nonpolar atoms are described by a 12-6 LJ potential. The LJ parameters are usually obtained from fitting to experimental vapor and liquid densities or vapor pressures. The Lorentz-Berthelot combining rules are used to describe the LJ interaction between different atom types8 σ RR þσββ ð1Þ σRβ ¼ 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi εRβ ¼ εRR 3 εββ ð1-kRβ Þ ð2Þ
^ 3alk , CO2^εalk , CO2 ψalk,CO2 3 m ^ alk m ^ CO2 σ nST,alk nST,CO2
¼
R
i¼1
¼
ncomp X ncomp X i¼1
j¼1
j¼1
xi xj
R¼1
β¼1
NR, i Nβ, j σ3Rβ εRβ
NR,alk Nβ,CO2 σ 3Rβ εRβ
ð4Þ
The parameters κRβ are optimized to accurately recover the value of the experimental property Ωexp in molecular simulation. In the remainder, we consider the simultaneous optimization of several cross interactions, denoted by κRβ. We define an objective function f(κRβ), which in the simplest case represents the leastsquared deviations of all observables, (Ωexp - Ωsim(κRβ))2. Within the force-field-based approach, the minimum of f(κRβ) cannot be found analytically; therefore, such a minimization has to be performed iteratively by changing the parameters κRβ.16-21 Here, we apply eq 5 to provide an estimate of how Ωsim(κRβ) depends on the parameters κRβ and obtain an analytic approximation for f(κRβ). In other words, we optimize the parameters κRβ to obtain dispersive cross-interaction energy parameters that result in the desired property Ωexp. The minimization of the objective function f(κRβ) with the PCP-SAFT equation of state only requires a few milliseconds on a desktop computer and enables the use of multidimensional optimizations (e.g., multiple temperatures or substances). As an example of this approach, the Henry coefficient Halk,CO2 of carbon dioxide in the homologous n-alkane series represents the observable Ωexp that we wish to predict by molecular simulations at varying temperature T. For details on these molecular simulations, we refer the reader to the Supporting Information. The solubility of carbon dioxide directly follows from the total pressure p with the Henry coefficient defined as ! yCO2 3 φVCO2 3 p Halk , CO2 ðTÞ lim ð6Þ xCO2 f0 xCO2
^ 3ij^εij ψij xi xj 3 m ^ jσ ^ im nST,i nST,j X X
β
We determine the value of ψij by requiring that the result of a molecular simulation for an observable Ωsim is best represented by the PCP-SAFT equation of state, as described in detail in the Supporting Information. The PCP-SAFT equation of state then becomes an explicit function of the force field parameters by regrouping eq 4 and combining with eq 2 as nST,alk nST,CO2 P P pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi NR,alk Nβ,CO2 σ3Rβ εRR 3 εββ ð1-kRβ Þ R β ^εalk ,CO2 ¼ ^ 3alk ,CO2 ψalk ,CO2 3 malk mCO2 σ ð5Þ
where we introduce a correction parameter κRβ to allow the cross-dispersive energy parameter εRβ to deviate from the geometric mean suggested by the simple Berthelot approximation. In the original TraPPE force field, κRβ = 0. When modeling the properties of mixtures with molecular simulations, one should adjust this correction parameter to reproduce experimental observables, such as bubble point pressures, enthalpies of vaporization, Henry coefficients, and so forth.6,14-21 In this way, one obtains a meaningful dispersive energy parameter εRβ for cross interactions. This is required, particularly, for modeling the interactions between polar and nonpolar atoms, as present in the n-alkane/carbon dioxide system.22,23 The PCP-SAFT equation of state applies the same methodology for cross-component interactions and offers a framework for the efficient optimization of the binary parameters κRβ. As polar interactions are included in PCP-SAFT, the appoach is valid for both polar and nonpolar systems. To approximate the value of an observable Ωsim obtained from molecular simulations with the PCP-SAFT equation of state, we need to establish a link between PCP-SAFT and force-field-based molecular simulation. The force-field-based approach can be related to the PCP-SAFT framework by equating their expressions for the Helmholtz energy. In the Supporting Information, we derive that for a system of ncomp components, this leads to nX comp n comp X
X X
ð3Þ
where the left-hand side of this equation is the leading dispersive term for the PCP-SAFT equation of state and the right-hand side equivalently represents the force-field-based model. The mole fraction of molecule type i is given by xi, and m̂i is its number of segments within the PCP-SAFT equation of state. Molecule type i consists of nST,i (united) atom types, and NR,i is the number of
where yCO2 is the carbon dioxide mole fraction in the vapor and φVCO2(T,p,yi) is the vapor-phase mixture fugacity coefficient of 394
dx.doi.org/10.1021/jz101688h |J. Phys. Chem. Lett. 2011, 2, 393–396
The Journal of Physical Chemistry Letters
LETTER
Figure 1. Molar Henry coefficient of carbon dioxide in butane. The solid diamonds represent the molecular simulation results with the TraPPE force field. The outcome from the simulations with the adjusted force field is given by the solid circles. Experimental data are given by the open triangles,24 and the lines are a guide for the eye. The errors in the computed Henry coefficients are smaller than the symbol size.
Figure 2. Molar Henry coefficient of carbon dioxide in decane. The solid diamonds represent the molecular simulation results with the TraPPE force field. The outcome from the simulations with the adjusted force field is given by the solid circles. Experimental data are given by the open triangles,24,25 and the lines are a guide for the eye. The errors in the computed Henry coefficients are smaller than the symbol size.
carbon dioxide. We assume that the binary parameters κRβ do not depend on the atom types of the carbon dioxide molecule, that is, κCH3,C = κCH3,O = κCH3,CO2 and κCH2,C = κCH2,O = κCH2,CO2. Starting from the TraPPE force field, the parameters κCH3,CO2 and κCH2,CO2 are then determined from the simultaneous leastsquares minimization of the PCP-SAFT results to experimental Henry coefficients for the n-alkanes series from C2 to C8 as well as decane, hexadecane, and eicosane according to eq 5, resulting in κCH3,CO2 = 0.053 and κCH2,CO2 = 0.101. These values are subsequently used in molecular simulations and are valid for the whole range of n-alkanes for 200 < T < 700 K. In principle, the minimum number of simulations needed to obtain values for κRβ equals the number of terms κRβ. In practice, it is convenient to fit these parameters to experimental data of several components. The results are illustrated for two representative systems, butaneCO2 (Figure 1) and decane-CO2 (Figure 2). These figures clearly show that the use of the binary parameters κCH3,CO2 and κCH2,CO2 leads to an excellent description of the experimental Henry coefficients over a wide temperature range, suggesting that the enthalpy of absorption is also correctly described.26 The original TraPPE force field (for which κCH3,CO2 = κCH2,CO2 = 0) clearly underestimates the experimental Henry coefficients. The average relative deviations between simulated and experimental Henry coefficients of butane and decane can be decreased from about 45% to approximately 3% by applying our new method. At high temperature, the temperature-dependent Henry coefficient passes through a maximum upon approaching the solvent critical temperature. This behavior is in agreement with other literature27 and is precisely captured by our corrected force field. In the Supporting Information, we show that by using the same values for κCH3,CO2 and κCH2,CO2, the solubility of carbon dioxide is accurately predicted for the full n-alkane series (C2-C20). This study proposes a simple method which requires a minimum number of simulations to provide cross-dispersion parameters that can quantitatively predict phase equilibrium properties. The corrected force field leads to accurate results for the solubility
of carbon dioxide in the n-alkane series C2-C20 and, thus, restores the transferability of the force field parameters for multicomponent systems. Most importantly, the method is not limited to alkane-carbon dioxide systems but is, in principle, valid for all systems that can be described by the PCP-SAFT equation of state. The method clearly shows the importance of combining molecular simulations and statistically based equations of state as efficient tools for the molecular design of solvents, for example, for finding more optimal solvents for carbon capture.
’ ASSOCIATED CONTENT
bS
Supporting Information. Detail description of how the PCP-SAFT equation of state and molecular simulations are related. We also provide details of our molecular simulations. This material is available free of charge via the Internet at http://pubs.acs.org.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ REFERENCES (1) Commission Staff Working Document SEC(2008) 85, Package of implementation measures for the EU’s objectives on climate change and renewable energy for 2020; 2008. (2) Metz, B.; Davidson, O. R.; Bosch, P. R.; Dave, R.; Meyer, L. A. Climate Change 2007: Mitigation. Contribution of Working Group Iii to the Fourth Assessment Report; Cambridge University Press: Cambridge, U. K., 2007. (3) Erdmenger, C.; Lehmann, H.; Mueschen, K.; Tambke, J. Klimaschutz in Deutschland: 40% Senkung der CO2 Emissionen bis 2020 Gegen€uber 1990; Umweltbundesamt Dessau: Dessau, Deutschland, 2007. (4) Zhang, L.; Siepmann, J. I. Direct Calculation of Henry’S Law Constants from Gibbs Ensemble Monte Carlo Simulations: Nitrogen, Oxygen, Carbon Dioxide and Methane in Ethanol. Theor. Chem. Acc. 2006, 115, 391–397. 395
dx.doi.org/10.1021/jz101688h |J. Phys. Chem. Lett. 2011, 2, 393–396
The Journal of Physical Chemistry Letters
LETTER
(5) Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A. H.; Ungerer, P. New Optimization Method for Intermolecular Potentials: Optimization of a New Anisotropic United Atoms Potential for Olefins: Prediction of Equilibrium Properties. J. Chem. Phys. 2003, 118, 3020– 3034. (6) Huang, Y. L.; Miroshnichenko, S.; Hasse, H.; Vrabec, J. Henry’s Law Constant from Molecular Simulation: A Systematic Study of 95 Systems. Int. J. Thermophys. 2009, 30, 1791–1810. (7) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244–1260. (8) Allen, M.; Tildesley, D. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (9) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569–2577. (10) Potoff, J. J.; Siepmann, J. I. Vapor-Liquid Equilibria of Mixtures Containing Alkanes, Carbon Dioxide, And Nitrogen. AIChE J. 2001, 47, 1676–1682. (11) Gross, J.; Sadowski, G. Application of the Perturbed-Chain SAFT Equation of State to Associating Systems. Ind. Eng. Chem. Res. 2002, 41, 5510–5515. (12) Gross, J. An Equation-of-State Contribution for Polar Components: Quadrupolar Molecules. AIChE J. 2005, 51, 2556–2568. (13) Gross, J.; Vrabec, J. An Equation-of-State Contribution for Polar Components: Dipolar Molecules. AIChE J. 2006, 52, 1194–1204. (14) Vrabec, J.; Stoll, J.; Hasse, H. Molecular Models of Unlike Interactions in Fluid Mixtures. Mol. Simul. 2005, 31, 215–221. (15) Vrabec, J.; Huang, Y. L.; Hasse, H. Molecular Models for 267 Binary Mixtures Validated by Vapor-Liquid Equilibria: A Systematic Approach. Fluid Phase Equilib. 2009, 279, 120–135. (16) Du, Q.; Yang, Z.; Yang, N. N.; Yang, X. N. Coarse-Grained Model for Perfluorocarbons and Phase Equilibrium Simulation of Perfluorocarbons/CO2 Mixtures. Ind. Eng. Chem. Res. 2010, 49, 8271–8278. (17) Baker, C. M.; Lopes, P. E. M.; Zhu, X.; Roux, B.; MacKerell, A. D. Accurate Calculation of Hydration Free Energies using PairSpecific Lennard-Jones Parameters in the CHARMM Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6, 1181–1198. (18) Schnabel, T.; Vrabec, J.; Hasse, H. Molecular Simulation Study of Hydrogen Bonding Mixtures and New Molecular Models for MonoAnd Dimethylamine. Fluid Phase Equilib. 2008, 263, 144–159. (19) Schnabel, T.; Vrabec, J.; Hasse, H. Unlike Lennard-Jones Parameters for Vapor-Liquid Equilibria. J. Mol. Liq. 2007, 135, 170–178. (20) Liu, A. P.; Beck, T. L. Vapor-Liquid Equilibria of Binary and Ternary Mixtures Containing Methane, Ethane, And Carbon Dioxide from Gibbs Ensemble Simulations. J. Phys. Chem. B 1998, 102, 7627– 7631. (21) Virnau, P.; M€uller, M.; MacDowell, L. G.; Binder, K. Phase Behavior of n-Alkanes in Supercritical Solution: A Monte Carlo Study. J. Chem. Phys. 2004, 121, 2169–2179. (22) Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z. Molecular Simulation of Phase Equilibria for Mixtures of Polar and Non-Polar Components. Mol. Phys. 1999, 97, 1073–1083. (23) Delhommelle, J.; Millie, P. Inadequacy of the LorentzBerthelot Combining Rules for Accurate Predictions of Equilibrium Properties by Molecular Simulation. Mol. Phys. 2001, 99, 619–625. (24) Yau, J. S.; Tsai, F. N. Solubilities of Carbon-Dioxide in NormalParaffins. Chem. Eng. J. Biochem. Eng. J. 1993, 52, 55–62. (25) Gmehling, J.; Rarey, J.; Menke, J. Dortmund Data Bank; Oldenburg, Germany, http://www.ddbst.com (2008). (26) Vlugt, T. J. H.; García-Perez, E.; Dubbeldam, D.; Ban, S.; Calero, S. Computing the Heat of Adsorption Using Molecular Simulations: The Effect of Strong Coulombic Interactions. J. Chem. Theory Comput. 2008, 4, 1107–1118. (27) Wilhelm, E. Low-Pressure Solubility of Gases in Liquids. In Measurement of the Thermodynamic Properties of Multiple Phases; Elsevier Science Bv: Amsterdam, The Netherlands, 2005; Vol. 7, pp 137-176. 396
dx.doi.org/10.1021/jz101688h |J. Phys. Chem. Lett. 2011, 2, 393–396
