6612
Ind. Eng. Chem. Res. 2007, 46, 6612-6629
Performance of COSMO-RS with Sigma Profiles from Different Model Chemistries Tiancheng Mu,†,‡ Ju1 rgen Rarey,† and Ju1 rgen Gmehling*,† Lehrstuhl fu¨r Technische Chemie, Institut fu¨r Reine und Angewandte Chemie, Carl Von Ossietzky UniVersita¨t Oldenburg, D-26111, Oldenburg, Germany, and Department of Chemistry, Renmin UniVersity of China, 100872, Beijing, People’s Republic of China
In this article, a comparison of the performance of COSMO-SAC and COSMO-RS(Ol), using σ profiles generated by different quantum mechanic programs (Gaussian 03, Turbomole and DMol3) and different model chemistries, with respect to activity coefficients at infinite dilution and vapor-liquid equilibrium data of binary systems, is presented. Density functional theory (BP and B3LYP) was used for the COSMO calculations. The basis sets used for the calculation were triple-ζ valence polarization (TZVP) and 6-311G(d,p). DMol3 results were taken from the publicly available VT 2005 databank. The predicted results were compared with the experimental data stored in the Dortmund Data Bank. The results showed that the σ profiles from B3LYP with 6-311G(d,p) basis set generated by Gaussian 03 lead to the best results for most binary systems; COSMORS(Ol) is superior to COSMO-SAC for predicting the thermodynamic properties of binary systems in most cases. Calculation results were analyzed, with respect to the types of components in the mixture. It can be expected that the performance of the models, when combined with B3LYP/6-311G(d,p), can be improved by adjusting the model parameters. 1. Introduction Knowledge of the thermodynamic properties of real mixtures (such as phase equilibria and excess properties) is a prerequisite for process design and product development. Such information is usually obtained from literature, experimental measurements, or prediction methods including corresponding states methods, equations of state, GE models, group contribution methods, quantitative structure-property relationship (QSPR) methods, Monte Carlo simulation, molecular dynamics methods, etc. The UNIFAC1 method, and its modified versions UNIFAC(Do)2,3 and UNIFAC(Ly),4 are the most reliable and widely used group contribution methods for engineers to predict the thermodynamic properties of mixtures. In group contribution methods, a pure compound or a mixture is considered as a mixture of predefined independent functional groups. The interaction energy of a system equals to the sum of group interaction energies, whereby the probability of certain groupgroup contacts is calculated via a Boltzmann expression. After the group interaction parameters are fitted to the experimental data of systems containing the same structural groups, the properties of any mixture composed of these structural groups can be calculated. Group contribution methods have been used to predict a variety of thermophysical properties such as vapor-liquid equilibria, liquid-liquid equilibria, solid-liquid equilibria, activity coefficients, excess enthalpies, azeotropic data, etc. They have also been applied to more-complicated systems, including polymer solutions, electrolyte solutions,5 etc. Besides the many advantages of group contribution methods, there are also some limitations. To predict the real mixture behavior, the interaction parameters between all groups in the * To whom correspondence should be addressed. Tel.: ++49-4417983831. Fax: ++49-441-7983330. E-mail address: gmehling@ tech.chem.uni-oldenburg.de. † Lehrstuhl fu ¨ r Technische Chemie, Institut fu¨r Reine und Angewandte Chemie, Carl Von Ossietzky Universita¨t Oldenburg. ‡ Department of Chemistry, Renmin University of China.
mixture must be known. Group contribution methods are not applicable to mixtures that contain new functional groups. In addition, interactions between groups inside one molecule and the position of the groups in the molecule are not considered; this leads to difficulties in regard to differentiating between isomers. An alternative and novel approach is to use solvationthermodynamics methods, utilizing surface charge densities calculated via the conductor-like screening model (COSMO), to characterize molecular interactions. The COSMO model is a dielectric continuum solvation model6,7 in which the dielectric constant of the continuum is set to infinity (i.e., an ideal conductor).8 This has the advantage that all shielding charges are located at the interface between the molecule and the continuum solvent has no electric field; hence, no charge distribution can exist inside the conductor. The solvation energy can be calculated from the surface shielding charges without the need to integrate over a certain distance into the solvent continuum. Based on this model, a conductor-like screening model for real solvents (COSMO-RS) was later developed by Klamt and co-workers.9-11 COSMO-RS also uses the local composition concept and allows the calculation of the chemical potential of the pure components and the components in mixtures. Instead of the potential of structural groups in the mixture, in COSMO-RS, the potential of a surface segment based on the shielding charge density of this segment and the other segments in the mixture is calculated. It is assumed that any segment can come into contact with every other segment. The calculation is based on the shielding charge density distribution of the surface elements of the molecules (σ profiles). However, when deriving the σ profiles, information on the mutual position of the segments, with respect to each other, is lost. To perform the COSMO-RS calculation, a small number of adjustable and universal parameters and the molecule-specific surface charge distribution are required. Some of these universal parameters can be obtained by measurement and others by regression to experimental data.10 After these parameters have
10.1021/ie0702126 CCC: $37.00 © 2007 American Chemical Society Published on Web 09/05/2007
Ind. Eng. Chem. Res., Vol. 46, No. 20, 2007 6613
been obtained, no further experimental data are required. They can even be used to calculate experimentally unavailable behavior, such as that of reactive intermediates or molecules in transition states. This, on the other hand, is also a disadvantage, because it is not easily possible to “train” the method by adjusting a larger number of parameters to the large amount of readily available information. Because no experimental data are required, a growing number of scientists and engineers show their interest to apply or further develop the COSMO-RS method. COSMO-RS can be used for the prediction of vapor-liquid equilibria (VLE), liquid-liquid equilibria (LLE), solid-liquid equilibria (SLE), vapor pressures of pure compounds and mixtures, partition coefficients, heats of vaporization, activity coefficients, solubilities, excess Gibbs free energies, and excess enthalpies, etc. COSMO-RS has also been extended to ionic liquids,12,13 polymers, surfactant micelles, biomembranes, and so on.14 The required COSMO calculation has been implemented in many quantum chemistry software packages such as Gaussian,15,16 Turbomole,17 MOPAC,18 DMol3,19 and GAMESS,20 etc. Two different further versions of COSMO-RS were also developed recently: COSMO-SAC (segment activity coefficient)21,22 and COSMO-RS(Ol)23. After the universal parameters are available, the σ profiles, surfaces, and volumes of the molecules involved in the system are the only properties required to perform the COSMO-RS calculation. σ profiles can be generated from different model chemistries implemented in various quantum chemistry software packages. In this paper, COSMO calculations based on Turbomole24-26 BP27,28/ TZVP,29 Gaussian 03 BP/TZVP and Gaussian 03 B3LYP/6-311G(d,p) were performed to obtain the charge density distribution on the molecular surface of compounds stored in the Dortmund Data Bank (DDB).30 To keep the results comparable, the same conformer was used in the different model chemistries. Activity coefficients at infinite dilution and VLE data were calculated using both the COSMO-SAC and COSMORS(Ol) methods, and these results were compared to experimental data stored in DDB. In addition, a comparison to the results of calculations using the σ profiles from the public VT 200531 databank (calculated with DMol3 using BP/DNP) is given, whereby the VT 2005 calculations are often for conformers that are different than the calculations performed within this work or supplied by DDBST, Ltd. [Here, the term “DNP” represents double numerical basis with polarization functions.] The detailed procedure of the calculation of σ profiles database VT 2005 has been explained by Mullins.31 A comparison to the COSMO-RS-version of Klamt and coworkers is not included here, because the authors are not able to exactly reproduce the calculation procedure used in the COSMOtherm software package. This version of COSMO-RS was only partly published in the open literature, and several shortcomings of the method discussed here may have been addressed by the numerous modifications in the commercial code. Different model chemistries usually lead to different results, and the COSMO-RS parameters should be adjusted to the model chemistry used: “COSMO-RS parameters are not specific of functional groups or molecule types. The parameters have to be adjusted for the QM-COSMO method that is used as a basis for the COSMO-RS calculations only.” (COSMOtherm Manual, C12 0103) In addition, for the calculation of σ profiles, a construction of the molecular-shaped cavity in the solvent continuum is required and surface segments need to be distributed over the
cavity surface. This procedure, as well as the required outlying charge correction, is implemented in different ways in the different quantum mechanical software packages, resulting in different calculation results, even in the case of identical atomic coordinates, model chemistries, and basis sets. To overcome this problem, COSMOlogic supplies specific parametrizations for the different quantum mechanical programs, as well as a conversion for calculations from Gaussian 03. In this work, we used the standard parametrizations of the different packages to determine the influence of the different implementations on the results of a COSMO-RS calculation, in comparison to the experimental findings. In addition, it was of interest whether a more modern model chemistry (B3LYP/6-311G(d,p)) would provide better results than BP/TZVP, which was used for the parametrization of both models. The B3LYP method yields more-accurate results for many properties, because of the use of a hybrid functional. The downside is that hybrid functionals cannot be combined with the RI (resolution of identity) method, which reduces memory requirements and significantly accelerates the calculation. 2. Theory and Model The concept for the description of the molecular interactions is based on the physical view of pairwise interacting surface segments. The difference between the screening charge densities σ and σ′ of a contact pair is a measure for the misfit of the screening charge density on both segments in the real system, compared to the situation inside the ideal conductor:
(R′2 )(σ + σ′)
Emisfit(σ,σ′) ) aeff
2
(1)
where aeff is the effective contact surface area and R′ is an energy factor.9,10 The surface charge density distribution of the molecules required for a COSMO-RS calculation is generated by the quantum mechanical calculation. An ensemble averaging procedure then is applied over a region of radius raν, using the following averaging algorithm:10
( ) ( ( ) (
∑µ σ*µ σν )
∑µ
rµ2raν2
rµ2 + raν2 rµ2raν2
rµ2 + raν2
exp -
exp -
dµν2 rµ2 + raν2 dµν2
rµ2 + raν2
)
)
(2)
where dµν is the distance between the segments µ and ν, rµ is the mean radius of segment µ calculated from the area of segment, and rRν is an adjustable parameter. The charge density distribution (σ profile, pi(σ)) is calculated from the averaged surface charge densities. The σ profile for a molecule i is defined as
pi(σ) )
Ai(σ) Ai
(3)
where Ai(σ) represents the total surface area of all the segments with the charge density σ and Ai is the total cavity surface area. Instead of a continuous function pi(σ), a discrete list of intervals pi(σn) is used (50 intervals in the range from -0.025 e/Å2 to +0.025 e/Å2 for the case of COSMO-SAC and 80
6614
Ind. Eng. Chem. Res., Vol. 46, No. 20, 2007
Figure 1. Comparison of σ profiles of acetone obtained from different model chemistries.
intervals in the range from -0.040 e/Å2 to 0.040 e/Å2 for the case of COSMO-RS(Ol)). In the case of COSMO-RS(Ol), a further averaging is applied to obtain the final σ profile:
pi(σn)COSMO-RS(Ol) )
1 n+1
∑ pi(σn)
3 n-1
(4)
In mixtures, the probability ps(σ) of finding a segment with a screening charge density σ is obtained from the mole fraction xi weighted sums of the σ profiles pi(σ) of all pure components i in the system.
ps(σ) )
∑i xiAipi(σ) ∑i
(5) xiAi
The steps to calculate the real mixture behavior using COSMO-RS are listed below: (1) The combinatorial part is calculated from the molecular surface Ai and the total cavity volume Vi reported by the quantum chemical software. Different combinatorial expressions are used for COSMO-RS(Ol) and COSMO-SAC. (2) The residual part is calculated from the segment activity coefficients, which must be computed via the solution of a selfconsistent equation. The detailed procedure, equations, and parameters used for COSMO-SAC and COSMO-RS(Ol) are given in refs 21 and 23. 3. Computational Details The molecular structures of the compounds were directly taken from the chemical structure database (ChemDB) of the DDB. To start from a reliable three-dimensional (3D) geometry for the COSMO calculations, first of all, a molecular mechanics optimization using the MOPAC program package has been performed for every compound. The COSMO method, as implemented in the Gaussian 0332 and Turbomole program package, is used for the quantum chemical calculation, together with the default parametrization. The COSMO calculations have been conducted using the density functional theory (DFT) levels
of BVP86 and B3LYP, combined with the basis sets of TZVP and 6-311G(d, p),33,34 respectively. For the cavity construction, the atomic radii from Klamt’s research, as delivered with the Turbomole and Gaussian 03 program package, were used directly, because these radii are not sensitive to the quantum-mechanical method that is used.11 The COSMO radii recommended by Klamt et al.35 should be ∼1.17 (( 0.02) times larger than the van der Waals radii reported by Bondi.36 Both the COSMO-SAC and COSMO-RS(Ol) calculations were performed using results from all quantum mechanical calculations, except in the case of the DMol3 calculations (VT 2005). This latter case is exempted because the comparisons of results are later only given for a common set of data, where results were available from all calculation options; including VT 2005 would have significantly reduced the size of this set of data. At the same time, VT 2005 results are often for conformers different than the results of the other quantum chemical calculations. 4. Results and Discussion 4.1. σ Profiles. Based on the theory and computational details, the σ profiles of 2014 compounds were obtained with Turbomole BP/TZVP, 1392 compounds were obtained with Gaussian 03 BP/TZVP, and 1363 compounds were obtained with Gaussian 03 B3LYP/6-311G(d,p). Using VT 2005, 1432 profiles that were calculated with DMol3 BP/DNP were available. Hereafter, in this paper, the following abbreviations are used for the different model chemistries:
SAC TM BP SAC G03 B3LYP SAC G03 BP SAC VT 2005 OL TM BP OL G03 B3LYP OL G03 BP OL VT 2005
COSMO-SAC using Turbomole BP/TZVP COSMO-SAC using Gaussian 03 B3LYP/6-311G(d,p) COSMO-SAC using Gaussian 03 BP/TZVP COSMO-SAC using DMol3 BP/DNP COSMO-RS(Ol) using Turbomole BP/TZVP COSMO- RS(Ol) using Gaussian 03 B3LYP/6-311G(d,p) COSMO- RS(Ol) using Gaussian 03 BP/TZVP COSMO-RS(Ol) using DMol3 BP/DNP
The σ profiles of acetone and n-hexane obtained from different calculations are shown in Figures 1 and 2, respectively. For the calculation, the COSMO-RS(Ol) charge density averaging
Ind. Eng. Chem. Res., Vol. 46, No. 20, 2007 6615
Figure 2. Comparison of σ profiles of hexane obtained from different model chemistries.
procedure was used in all cases. The charge distributions calculated with different models obviously are not identical. To analyze the effect of these variations on the calculated mixture behavior, activity coefficients at infinite dilution and VLE data were calculated and compared to the experimental data stored in the DDB. 4.2. Results for Activity Coefficients at Infinite Dilution. The activity coefficients at infinite dilution (γ∞) are of great practical importance for the simulation of distillation processes, environmental protection, etc. At the same time, they usually represent the highest deviation from ideal behavior. Various techniques for the measurement of γ∞ are available, and more than 46 000 values have been published in the literature. The modified UNIFAC(Do) procedure is commonly used to predict these data, because γ∞ data were used to regress the group interaction parameters. COSMO-RS presents an alternative in case the required interaction parameters are not available.37 4.2.1. Overall Root-Mean-Square Deviations (RMSD) and Relative Average Deviations in γ∞. The γ∞ coefficients of binary systems were calculated using COSMO-SAC and COSMO-RS(Ol), with σ profiles generated using the different software and model chemistry options. The root-mean-square deviations (RMSD) of the calculated results, with respect to the experimental data, were calculated using the following equation:
RMSD )
x∑ 1 n
(ln(γ∞calc) - ln(γ∞exp))2
(6)
n
where n is the number of data points and ln(γ∞calc) and ln(γ∞exp) are the respective logarithms of the calculated and experimental activity coefficients at infinite dilution. The RMSD values are depicted graphically in Figure 3. More than 46 000 experimental γ∞ values are stored in the DDB (Version 2006). The published γ∞ data cover the range from 0.02 to 1010. Usually, the large values of γ∞ are afflicted with larger experimental errors; the data measured by different authors sometimes differ by more than 1 order of magnitude. Therefore, for the model comparison, the data base was limited to γ∞ values of
