COSMO-RS Methods

In principle, all electronic and molecular interactions are described by quantum mechanics, so you may wonder why we have not considered computing mixture properties from this fundamental approach. In practice, two considerations limit the feasibility of this approach. First, quantum mechanical computations tend to be time consuming. Precise computations can require days for a single molecule the size of naphthalene and the computation time increases as N⁷, where N is the number of atoms in the molecule. Second, the intermolecular interaction energy, which affects the mixing properties, is at least an order of magnitude smaller than the intramolecular energy. Therefore, high precision is required to compute the intermolecular interactions directly. Circumventing these limitations has been the focus of the “COSMO-RS” approach.

COSMO-RS is an abbreviation for “COnductor-like Screening MOdel for Real Solvents.” It refers to a method of performing quantum mechanical calculations as if the simulated molecule were in a conductive bath rather than a vacuum. The method was developed originally by Klamt and coworkers¹⁵ as an extension of previous work on a continuum solvation model (CSM).¹⁶ The implementation of Klamt and coworkers is marketed commercially as COSMOtherm and updated continually. A free educational version with graphical user interface is available that includes capability for about 350 compounds. The implementation of Klamt and coworkers is currently based on the TURBOMOLE package for quantum mechanical simulations, and a few other packages also provide consistent results. Later, Lin and Sandler developed an implementation based on the DMOL³ simulation package,¹⁷ which is included in the Accelrys Materials Studio. We refer to this implementation as COSMO-RS/SAC. It is available as a free download including capability for about 1500 compounds from a web site maintained by Y.A. Liu at Va. Tech.¹⁸ The method includes several empirical parameters that have been characterized by fitting experimental data. The specific parameters depend on the quantum simulation method. In the example below, we have applied the SAC parameters.

The key idea of the COSMO-RS approach has been to focus on the polarization of the surface surrounding a molecule. The significance of the surface polarization can be likened to the acidity and basicity of the SSCED model. If the surrounding surface is positively charged, then the molecular surface must be basic; if the surrounding surface is negatively charged, then the molecular surface must be acidic. If we imagine coloring these acidic interactions as blue and the basic interactions as red, we could represent the molecular surfaces in the manner of the pictures on the cover of the textbook. The overall surface energy, including positive, negative, and neutral influences, can be correlated with activity coefficients and other properties like the heat of vaporization.

To calculate the surface polarization an approximate quantum mechanical method can be applied, known as density functional theory (DFT). DFT computations typically require only a few hours for a molecule like naphthalene, and the increase in computation time scales as N³. If we integrate over the interactions between the surfaces of molecules, we can imagine how results similar to SSCED could be achieved. An additional advantage of COSMO-RS is that local composition effects are implicit in the integration of the local polarization interactions over all orientations between the two molecules. Furthermore, acidity and basicity are not limited to a single characteristic value per molecule, but are characterized by a range of polarizations over the entire surface of each molecule.

To implement the method, the molecular surface charges are calculated using DFT. The observed polarization is referred to as a σ-profile, where σ refers to the surface charge density (Coulombs/Ų). Typical σ-profiles are illustrated in Fig. 13.3. The p(σ) represents the amount of area per σ-interval (Ų/(Coulombs/Ų)) plotted against the surface charge density (Coulombs/Ų). In other words, p(σ) is proportional to a count of how many segments possess a given amount of surface polarization. It is analogous to the number of occurrences of a group in UNIFAC. The curve is normalized such that an integral of the σ-profile gives the total surface area of the molecule. The effective area per segment is divided out near the end of the calculation when computing the activity coefficient. The area under the curve in a particular region shows the total area of the molecule with the charge density. Fig. 13.3(a) shows how ethanol is a smaller molecule than octanol based on the total area under the curve.

Figure 13.3. Samples of σ-profiles for application to the COSMO-RS method, A_effn_i(σ). Dashed vertical lines show the threshold values for hydrogen bonding.

Since the charge density, σ, ranges over negative and positive values, all values less than –0.0084 are considered to contribute to acidity. A similar consideration applies to the β contribution except that all σ > 0.0084 contribute to basicity. In Fig. 13.3(a), note that the two alcohols have extremely similar contributions of total area for both acidity and basicity. Fig. 13.3(b) shows that water is a very small molecule (small area under the curve) with extensive hydrogen bonding (both acidic and basic outside the dispersion bounds) while chloroform is a relatively large molecule with strong acidity and no basicity. All of these behaviors make sense qualitatively, so what remains is the translation into a quantitative method.

Computing activity coefficients for a binary solution is similar to the UNIFAC method if you can imagine transforming the earlier summation over groups to a summation over discretized polarization segments of the σ-profile. Note that use of the term “segment” does not refer to a geometric segment, but to an “interval” of surface polarization. A particular amount of surface polarization may occur at various places over the surface of the molecule, but all would be added together to get p(σ). It may be helpful to think of this quantity in mathematical terms rather than as a physical entity. Typically, the σ-profile is discretized into 51 values ranging from –0.025 to 0.025 (Coulombs/Ų).¹⁹ Therefore, we can refer to σ_k where k = [1,51]. To simplify the computations, the integration ∫p(σ)dσ ~ ΔσΣp(σ_k) is represented as a sum, Σp(k). The larger size of a particular molecule is reflected in the Σp(k) in a comparable manner to the molecular volume in the Scatchard-Hildebrand or SSCED models. Discretized segments on molecule i interact with segments on all molecules, including the i^th molecule itself. In a mixture, p_i(k) designates p(k) for the i^th component. The activity coefficient contribution for a segment k(Γ_k) in a solution of polarization segments is analogous to Eqn. 13.57, but with a significantly different functional form.

The surface area fraction of polarization segments in a binary mixture of molecules type 1 + 2 is

where a_jk is defined below and the molecular surface area is

A practical difference between Eqns. 13.57 and 13.61 is that the activity coefficient contribution of a given σ-interval depends on itself through the summation over all interactions in Eqn. 13.61. Conceptually, the influence of a segment on its neighbors alters its own behavior as those neighbors respond to the local activity. This necessitates iteration, initiated with Γ_k = 1 for all k.

A fundamental difference between Eqns. 13.57 and 13.61 is that the formulas are derived entirely differently. The matrix of a_jk is related to balancing electrostatic charges, as given by:

where α(j,k) and β(j,k) characterize the hydrogen bonding between the j^th and k^th segments, similar to the α and β parameters of the SSCED model. We assume that hydrogen bonding occurs regardless of which σ value (j or k) surpasses the threshold value because the energetic reward is sufficient for them to find each other regardless of where the segments are. Mathematically, this becomes,

Similar to UNIFAC, this approach leads to a nonunity value for the activity coefficient of a pure fluid. So,

where  is computed using the same formula as for UNIFAC or UNIQUAC; A_eff = 7.5 is the normalization for the area, and Γ_k⁽ⁱ⁾ is the segment activity in pure component i, computed by applying Eqn. 13.62 and so forth, with the appropriate composition.

Although we can conceive of COSMO-RS as being similar to UNIFAC, it is really much more general. The UNIFAC method requires experimental data to characterize the a_ik matrix. COSMORS, on the other hand, computes this matrix based on σ-profiles that have been computed with no experimental data except the empirical constants in Eqn. 13.65. Usually, more precise results are obtained with UNIFAC if all the groups have been characterized, but the COSMO-RS approach can be used to supplement the UNIFAC method when no experimental data exist.

---

Example 13.7. Calculation of activity coefficients using COSMO-RS/SAC

σ-profiles for methanol and acetone are listed below. (a) Use these to compute the activity coefficient at x_M = 0.425 and 55°C assuming the SAC values of the COSMO-RS parameters as given by Lin and Sandler. (b) Compute the activity coefficients over the entire range of compositions and compare to the fit of UNIQUAC to experimental data and the SSCED model.

Table of σ-profiles. Note that σ(k) = [-0.025,0.025], a total of 51 values. Cells omitted if all zeros.

Solution

a. Applying Eqn. 13.64 gives q₁ = 67.9 and q₂ = 102.6. The mixture’s segmental area fraction, Θ(k), is zero for k = 1 to 9. Θ(10)=(0.425(0.598) + (0.575)0)/(0.425(67.9) + (0.575)102.6) = 2.89(10^-3) and so forth for Θ(11) to Θ(42).

For j = k = 1, σ(j) = σ(k) = –0.025. So, max(σ(j), σ(k)) = –0.025, and –0.025 – 0.0084 = – 0.0334, but max(0, –0.0334) = 0.

So α(1, 1) = 0, and even though β(1, 1) = –0.0166, the product α(1, 1)β(1, 1) =0 and hydrogen bonding contributes zero for i = k = 1. For the dispersion term, however,
a_1,1 = [8233(σ(1) + σ(1))²]/0.001987 = 10359 and Ψ(1,1) = exp(–10359/328.15)=1.951(10^–14).

The first term for which hydrogen bonding is nonzero is i = 1, k = 35. Then, α(1, 35) = max(σ(j), σ(k)) – 0.0084 = 0.0006 and β(1, 35) = –0.0166.

The total is a_1,35 = [8233(–0.025 + 0.009)² + 85580(0.0006)(–0.0166)]/0.001987= –428.98 and Ψ(1, 35) = exp(428.98/328.15) = 3.696. Note how p(k) contributes to Θ but not Ψ.

These two cases (i = k = 1) and (i = 1, k = 35) illustrate important behaviors. Briefly, segment pairs involving the disperse interactions between two like-charged segments are disfavored, as indicated by the small value of Ψ. Even when the polarizations are exactly opposite and sum to zero, the best that can happen for the dispersion term is Ψ = 1. All other values of polarization lead to Ψ<1 for the dispersion term (because the contribution is squared), reflecting an unfavorable interaction. This is vaguely reminiscent of the (Δδ’)² term that applies to dispersion interactions in the SSCED model. The second case shows that acid-base interactions are opposite in sign to the dispersion interactions, and generally are sufficient to provide Ψ > 1, indicating favorable interactions. Again, this is reminiscent of hydrogen bonding in the SSCED model.

Completing the solution for part (a) involves summing all the terms, both for the mixture and for the pure fluid. After the first iteration of Eqn. 13.61 for Γ_k, the first 10 values are roughly 0.5. Other sample values are Γ₂₀ = 1.561, Γ₃₀ = 1.203, Γ₄₀ = 0.736, and Γ₅₀ = 0.500. After the last iteration, these values become Γ₁ = 7.864(10^–5), Γ₁₀ = 0.0279, Γ₂₀ = 1.478, Γ₃₀ = 1.107, Γ₄₀ = 2.522, Γ₅₀ = 3.312(10^–3).

Similarly, after the first iteration for pure fluids: Γ_k⁽¹⁾ = 0.5 for k = 1 to 10, …, Γ₂₀⁽¹⁾ = 1.559, Γ₃₀⁽¹⁾ = 1.267, … Γ₁⁽²⁾ = 0.501, Γ₁₀⁽²⁾ = 0.563, Γ₂₀⁽²⁾ = 1.562, and Γ₃₀⁽²⁾ = 1.176.

After the final iteration, Γ₁⁽¹⁾ = 7.05(10^–5), Γ₁₀⁽¹⁾ = 2.931(10^–2), Γ₂₀⁽¹⁾ = 1.546, Γ₃₀⁽¹⁾ = 1.204, … Γ₁⁽²⁾ = 3.322(10^–4), Γ₁₀⁽²⁾ = 3.748(10^–3), Γ₂₀⁽²⁾ = 1.453, Γ₃₀⁽²⁾ = 1.063, …

Summing the Γ_k – Γ_k⁽ⁱ⁾ gives lnγ₁ = 0.0957 – 0.0353 = 0.0605 and lnγ₂ = 0.0746 – 0.0145 = 0.0601, where the negative terms are the Guggenheim-Staverman contribution using q_i and r₁ = 48.8 and r₂ = 86.4. The values of r₁ and r₂ come from the (VT) database of Mullins et al.

b. The σ-profiles for methanol and acetone are shown below, along with a comparison to the activity coefficients from the SSCED model and the UNIQUAC model.^a

Regarding the σ-profiles, it is apparent from the high positive polarization (proton acceptor) of the acetone that hydrogen bonding should play a significant role. The right figure shows the comparison of the COSMO-RS model with a UNIQUAC fit (considered the benchmark in this case). The COSMO-RS/SAC model and the SSCED model are almost scaled mirror images, and neither is precisely correct.^b The COSMO-RS/SAC model underestimates the nonideality of this mixture.

---

a. The UNIQUAC model (right figure) was fit to the data of Marinichev A.N., Susarev M.P. 1965. Zhur. Prikl. Khim, 38(2):371.

b. These components are available in the educational version of COSMOtherm. Estimates are more accurate with COSMOtherm, but the mirror-image effect is consistent.