Statistical Thermodynamic Models: COSMO-RS, COSMO-SAC

The next sections will give an overview of different thermodynamic models ranging from implicit solvation models to Non-random Two-liquid (NRTL) and group contribution (GC) models. Many of these models found on common theory and can even be combined to obtain fine-tuned parameter sets.

Currently the Cebule engine supports the following models:

• COSMO-RS & COSMO-SAC (already supported in the SDK) and various extensions: COSMO-RS + dispersion, COSMO-SAC-dsp, eCOSMO-RS/-SAC

• Non-random Two-liquid (NRTL, eNRTL) and Wilson models (upcoming integration in SDK)

• Group contribution (GC) models including modified UNIFAC (Dortmund; supported by SDK), UNIFAC (Lyngby; upcoming integration in SDK), UNIFAC 2.0 (Kaiserslautern; upcoming integration in SDK)

• Hybrid models: UNIFAC-VISCO (supported by SDK), UNIFAC-VISCO-IL (supported by SDK) and COSMO-NRTL (upcoming)

On this page we first present how the software development kit (SDK) can be utilized to retrieve sigma profiles and sigma-moments from the COSMO-RS and COSMO-SAC models, different parameterisations and model extensions. The theoretical theory of the models is covered after the SDK description.

Create a COSMO-RS/-SAC task with the Cebule SDK

Cebule SDK TaskType: SIGMA

• Calculates sigma profiles from COSMO surface charge densities
• Requires: A completed COSMO task (must be connected via connected_task_id)
• Inputs:
• cosmo_method: str - Method for calculating sigma profiles, from [cosmo-rs, cosmo-sac]
• For COSMO-SAC only:
  • Optional: num_profiles: int - Number of profiles to compute (1 or 3, default: 3 and 1 only works in certain molecules)
  • Optional: averaging: str - Averaging method (default: "Mullins"), from [Hsieh, Mullins]
• Cebule max_processors: Not used
• Output: Dictionary containing sigma profiles and related information:
• For COSMO-RS:
  • Charge bins and charge densities across the molecular surface
  • Sigma moments list (6 sigma moments)
  • M₀: Molecular surface area of the compound
  • M₁: Negative of total charge of the compound
  • M₂: Polarity of the compound
  • M₃: Asymmetry of the sigma profile
  • M_4-6: Higher-order moments (without well-established physical interpretations)
  • Hydrogen bond donor and acceptor sigma moments lists (6 sigma moments in each)
  • M_hbacc: Quantifies the molecule's hydrogen bond acceptor strength
  • M_hbdon: Quantifies the molecule's hydrogen bond donor strength
• For COSMO-SAC:
  • Three separate sigma profiles for distinct surface types:
  • NHB: Non-hydrogen bonding surface segments. Generally correspond to hydrophobic or weakly polar regions of the molecule.
  • OH: Hydroxyl group surface segments. These areas have distinctive charge density distributions that reflect their strong hydrogen bond donor/acceptor capabilities.
  • OT: Other hydrogen bonding surface segments. This includes N, F, and O (outside of OH) atoms

Example:

# First create a COSMO calculation
cosmo_task = session.cebule.create_task("COSMO Ethanol",
                                     TaskType.COSMO,
                                     cosmo_input,
                                     nprocs=16)

# Then calculate sigma profiles using the COSMO results
sigma_task = session.cebule.create_task("Sigma Profile Ethanol",
                                     TaskType.SIGMA,
                                     connected_task_id=cosmo_task.id,
                                     cosmo_method="cosmo-sac")

# For COSMO-SAC with custom parameters
sigma_task_custom = session.cebule.create_task("Sigma Profile Custom",
                                           TaskType.SIGMA,
                                           connected_task_id=cosmo_task.id,
                                           cosmo_method="cosmo-sac",
                                           num_profiles=3,
                                           averaging="Hsieh")

# For COSMO-RS
sigma_task_rs = session.cebule.create_task("Sigma Profile COSMO-RS",
                                        TaskType.SIGMA,
                                        connected_task_id=cosmo_task.id,
                                        cosmo_method="cosmo-rs")

Cebule SDK TaskType: SOLUBILITY

• Calculates solubility of a solute in solvent system using COSMO-SAC or COSMO-RS sigma profiles
• Requires: Completed SIGMA tasks for both solute and solvent(s) (must be connected via connected_task_id). Pass all connected task IDs in a list, with the solute task first.
• Inputs:
• temperature: float - Temperature in Kelvin
• melting_point: float - Melting point of the solute in Kelvin
• enthalpy_melting: float - Enthalpy of fusion in J/mol
• change_heat_capacity_melting: float - Heat capacity change upon melting in J/(molK). Note*: Only used for COSMO-SAC sigma profiles; ignored for COSMO-RS sigma profiles
• sol_init: float - Initial guess for solubility value
• solv_composition: List[float] - Molar fractions of solvents in a mixed solvent system. Should be in the same solvent component order as the connected task IDs (exclusing the first connected task ID which is for the solute).
• Cebule max_processors: Not used
• Output: Dictionary containing solubility results:
• solubility: float- Predicted solubility value (mole fraction)
• Optional: uncertainty: float- Uncertainty on the predicted solubility value, for COSMO-SAC

Example:

# First create SIGMA tasks for solute and solvents
# Assuming we already have COSMO tasks for each molecule
sigma_solute = session.cebule.create_task("Sigma Profile Solute",
                                       TaskType.SIGMA,
                                       connected_task_id=cosmo_solute_task.id,
                                       cosmo_method="cosmo-sac") # Can also use "cosmo-rs"

sigma_solvent1 = session.cebule.create_task("Sigma Profile Solvent1",
                                         TaskType.SIGMA,
                                         connected_task_id=cosmo_solvent1_task.id,
                                         cosmo_method="cosmo-sac")

sigma_solvent2 = session.cebule.create_task("Sigma Profile Solvent2",
                                         TaskType.SIGMA,
                                         connected_task_id=cosmo_solvent2_task.id,
                                         cosmo_method="cosmo-sac")

# Calculate solubility (first sigma task is assumed to be the solute)
solubility_task = session.cebule.create_task("Solubility Calculation",
                                          TaskType.SOLUBILITY,
                                          connected_task_id=[sigma_solute.id, 
                                                           sigma_solvent1.id, 
                                                           sigma_solvent2.id],
                                          temperature=298.15,
                                          melting_point=423.15,
                                          enthalpy_melting=25000.0,
                                          change_heat_capacity_melting=35.0, # Only specify this for COSMO-SAC
                                          sol_init=0.01,
                                          solv_composition=[0.7, 0.3])  # 70% solvent1, 30% solvent2

The SOLUBILITY task leverages COSMO-SAC/RS sigma profiles to predict the solubility of compounds in pure solvents or solvent mixtures.

1. The first connected SIGMA task is treated as the solute
2. The remaining connected SIGMA tasks are treated as solvents
3. The solv_composition parameter must match the number of solvents
4. All sigma profiles must be either all COSMO-SAC profiles, or all COSMO-RS profiles

This approach provides accurate solubility predictions for drug-like compounds and other organic molecules in common solvents, which is valuable for formulation development, crystallization process design, and other pharmaceutical and chemical applications.

COnductor like Screening MOdel for Real Solvents (COSMO-RS)

The COSMO-RS model invented by Andreas Klamt is founded on statistical thermodynamics and utilises information from quantum chemical calculations. In the COSMO-RS model, no molecular fragments/groups or component specific parameters are necessary to be fitted to experimental data as is the case for the UNIFAC model. Although the COSMO-RS model parameter set has also been fitted to experimental data, theoeretically any molecular structure can be calculated with COSMO-RS whereas with the UNIFAC model a compound can not be predicted if the experimental data for fitting the model has not included a specific molecular fragment of the molecule of interest in the data set.

The value of COSMO-RS to the pharma, life science and chemical industry is its predictive and computational performance by performing an ab-initio density functional theory (DFT) level calculation for each compound of interest. The generated data from the implicit solvation model COSMO (cosmo data) and from the post-process algorithm in form of the COSMO-RS model (sigma data) allows one to store a lookup data table for each compound to then efficiently calculate thermodynamic properties such as solubility and phase equilibria. Further, the generated descriptor data in form of sigma-profiles or sigma-moments allow to be utilized for classification and regression models (QSPR, machine learning etc.) and can be combined with other descriptors from for example molecular dynamics (MD) calculations to obtain highly predictive QSPR/ML models for viscosity and other important thermophysical properties.

The fact that the COSMO-RS model relies on stored cosmo data and the thermodynamic property prediction step relies on a lookup table of sigma profiles allows to calculate temperature dependent properties as a sub-routine in high performance iterative algorithms such as computational fluid dynamics (CFD) models or process simulation software. These kind of models/softwares need to be able to receive property estimates in a timely manner (within milliseconds to seconds). An iterative method which would include a DFT calculation at each step would not be feasible for such kind of models, algorithms and softwares.

COSMO-RS is based on the assumption that the interaction between molecules can be represented as an ensemble of interacting surface segments :

\[ \mu_S(\sigma_i) = -kT ln (\int \sigma_j p_s(\sigma_j) exp[(-E_{ij} + \mu_s(\sigma_j))/kT])\]

The quantum chemical calculation used to obtain these surfaces are performed in the conductor state, a state in which a molecule makes up a cavity surrounded by a liquid medium which is assumed to be perfectly conducting (\(\epsilon \rightarrow \infty\)).

The value of the surface charge density from COSMO (CPCM) calculations is averaged to reduce numerical noise: $$ \sigma_{\alpha} = \frac{\sum_{\beta}\sigma_{\beta}^*\frac{r_{\beta}^2 r_{av}^2}{r_{\beta}^2 + r_{av}^2}exp\left(\frac{-d_{\alpha \beta}^2 }{r_{\beta}^2 + r_{av}^2}\right)}{\sum_{\beta}\frac{r_{\beta}^2 r_{av}^2}{r_{\beta}^2 + r_{av}^2}exp\left(\frac{-d_{\alpha \beta}^2 }{r_{\beta}^2 + r_{av}^2}\right) } $$

The surface segments are subsequently distributed onto a predefined grid of screening charge density values. As the \(\sigma\) values are situated between two gridpoints, they are projected onto both gridpoints, and the segment size is scaled with respect to the distance to each of the neighbouring grid points.

Assuming a fixed area of segment contact, the fraction of a segment type I in a mixture can be expressed as:

\[ X_{I} = \frac{\sum_i x_i A_I^i}{\sum_{i}\sum_{J} x_i A_J^i} \]

Where \(A_I^i\) is the area of segment type I on molecule i and \(x_i\) is the mole fraction of molecule i.

Apart from the screening charge density itself, a descriptor for the polarity of surrounding segments is added, namely the correlation screening charge density \(\sigma^{\perp}\).

\[ \sigma_{\alpha}^{\perp} = \sigma_{\alpha}^0 -0.816 \sigma_{\alpha} \]

To obtain \(\sigma_{\alpha}^0\) an additional averaging is performed with twice the averaging radius.

The contribution to the electrostatic misfit interaction energy to the interaction free energy, \(G^{int}_{IJ}\), is calculated as:

\[ E^{mf}_{IJ} = 0.5 a_{eff} \alpha ((\sigma_I + \sigma_J)^2 +f_{corr}(\sigma_I + \sigma_J)(\sigma^{\perp}_I + \sigma^{\perp}_J)) \]

Where \(\alpha\) is a misfit prefactor, \(a_{eff}\) is the area of segment contact, and \(f_{corr}\) is the correlation correction factor parameter that was adjusted based on the dielectric screening energies of a set of diverse molecules.

The hydrogen bond interaction energy (hydrogen bond contribution), \(G_{IJ}^{hb}\), is calculated as:

\[ G_{IJ}^{hb} = \begin{cases} c_{hb}(T)a_{eff}[min(0;f_{e_I}(\sigma_I+\sigma_{hb})) \cdot D_I \cdot max(0;f_{e_J}(\sigma_J-\sigma_{hb})) \cdot A_J] & \text{if $\sigma_I < \sigma_J$ }\\ c_{hb}(T)a_{eff}[ min(0;f_{e_J}(\sigma_J+\sigma_{hb})) \cdot D_J \cdot max(0;f_{e_I}(\sigma_I-\sigma_{hb})) \cdot A_I ] & \text{if $\sigma_I \ge \sigma_J$ } \end{cases} \]

Where \(\sigma_{hb}\) is the hydrogen bond threshold parameter and the \(f_{e_I}, f_{e_J}\) parameters give the opportunity for element specific scaling of the hydrogen bond strength. This parameter is set to unity by default. The min() function in the above equation represents the donor contribution and the max() function represents the acceptor contribution. The descriptors \(A_I, A_J\) and \(D_I, D_J\) allow to distinguish between segments from hydrogen atoms with negative screening charge densities (donors) and segments from potential hydrogen accepting atoms with positive screening charge densities (acceptors).

The descriptor definitions are:

\[ \begin{split} D_I &= \begin{cases} 1 & \text{if $\sigma_I < 0$ and $e_I$ in $(1,[100-150])$ }\\ 0 & \text{otherwise} \end{cases} \\ A_I &= \begin{cases} 1 & \text{if $\sigma_I \ge 0$ and $e_I$ in $(6,7,8,9,15,16,17,35,53)$ }\\ 0 & \text{otherwise} \end{cases} \end{split} \]

Where \(e_I\) is the element number of segment I.

The temperature dependence of the hydrogen bond interaction free energy parameter \(c_{hb}(T)\) is:

\[ c_{hb}(T) = c_{hb} \left[ max\left(0;1-c_{hb}^T + \frac{c_{hb}^T * 298.15}{T}\right) \right] \]

The interaction free energy \(G^{int}_{IJ}\) is restricted to the hydrogen bonding contribution and the misfit contribution. The complete interaction free energy is:

\[ G^{int}_{IJ} = E^{mf}_{IJ} + G^{hb}_{IJ} \]

The interaction parameter \(\tau_{IJ}\) is defined as:

\[ \tau_{IJ} = exp\left(\frac{-G^{int}_{IJ}}{RT}\right) \]

An important part of the COSMO-RS model is a non-linear set of equations named COSMOSPACE ¹. It relates the interaction free energies of the segment to an activity coefficient, \(\Gamma_I\), in a mixture of contacting segments:

\[ \Gamma_{I} = \frac{1}{\sum_{J}X_J\Gamma_J\tau_{IJ}} \]

The equation is solved iteratively until the result converges with a defined step size (openCOSMO-RS default: \(0.7\Delta \Gamma\)).

The logarithmic molecule activity coefficients with respect to the COSMO reference state can be calculated via:

\[ ln(\gamma_i^{res} (\textbf{x})) = \sum_I \frac{A^i_I}{a_{eff}} ln(\Gamma_I) \]

where \(\textbf{x}\) is the composition (mole fractions of the molecules) in the mixture. To retrieve the residual contribution of the molecular activity coefficient, the result must be converted to the pure component reference state by subtracting the molecular activity coefficient in the chosen reference state with the molecular activity coefficient in the COSMO reference state.

Thus, one obtains the residual molecular activity coefficient:

\[ ln(\gamma_i^{res}(\boldsymbol{x})) = ln(\gamma_i^{res}(\boldsymbol{x}))- ln(\gamma_i^{res}(\boldsymbol{x}^{ref})) \]

which can also be written as:

\[ ln \gamma_i^{res} = \int p^{(i)}(\sigma) |\mu_S(\sigma) - \mu_{S=(i)}(\sigma)| d\sigma \]

Three different options for the reference states can be set: Pure, reference mixture and cosmo. If for example cosmo is chosen no conversion to a different reference state is performed:

\[ ln(\gamma_i^{res}(\boldsymbol{x})) = ln(\gamma_i^{res}(\boldsymbol{x})) \]

The combinatorial term is calculated by using the Staverman-Guggenheim term ²:

\[ \frac{\phi }{\chi_i } = \frac{V_i}{\sum_j \chi_j V_j} \]

\[ \frac{\theta }{\chi_i } = \frac{A_i}{\sum_j \chi_j A_j} \]

\[ ln(\gamma_i^{comb}(\boldsymbol{x})) = ln\left(\frac{\phi}{\chi_i} \right) + 1 - \frac{\phi}{\chi_i} - \frac{z}{2} \frac{A_i}{A_{std}} \left( ln\left( \frac{\phi}{\theta}\right) +1 - \left( \frac{\phi}{\theta}\right) \right) \]

The combinatorial term is converted to the appropriate molecular reference state:

\[ ln(\gamma_i^{comb}(\boldsymbol{x})) = ln(\gamma_i^{comb}(\boldsymbol{x}))- ln(\gamma_i^{comb}(\boldsymbol{x}^{ref})) \]

with \(V_i\) and \(A_i\) being the volume and surface area (relative to a \(CH_2\) molecular fragment) of molecule i with mole fraction \(\chi_i\), z being the combinatorial number (number of close interacting molecules around the center molecule) here chosen as 10. \(A_i\) is the surface area of compound i and \(A_{std}\) is the normalized surface area (\(A_{std}\) is defined as q for COSMO-SAC(2010) by Hsieh, Sandler and Lin ³ and has been taken over by some COSMO-RS parameterizations).

The total activity coefficient is then calculated with:

\[ ln(\gamma_i(\boldsymbol{x})) =ln(\gamma_i^{comb}(\boldsymbol{x})) + ln(\gamma_i^{res}(\boldsymbol{x})) \]

In addition to calculating the activity coefficient of molecule i, the openCOSMO-RS implementation of the COSMO-RS model also allows one to calculate the contact probability between segments of type I and J. The contact probabilities are normalized with respect to segment type I where the resulting matrix is not symmetric:

\[ \frac{P_{IJ}}{X_I} = X_J\Gamma_I\Gamma_J\tau_{IJ} \]

The segment contact probability can be used to calculate the internal energy \(U\) of the system:

$$ U = \frac{N}{2} \sum_I \sum_J P_{IJ} E_{IJ} $$ where \(E_{IJ}\) is the interaction energy from a segment pair of type (I-J) and N is the number of segments in the mixture. The average interaction energies of the surface segments of molecules of type i in the mixture can be calculated from:

\[ \epsilon_i = \sum_I \sum_J N^i_I \frac{P_{IJ}}{X_I} E_{IJ} \]

where \(N_I^i = \frac{A_I^i}{a_{eff}}\)

Benchmarking different parameterisations of COSMO-RS/-SAC model versions and extensions

The aim of this detailed documentation of the COSMO-RS and COSMO-SAC models, their extensions and parameterisations, is to also lay the foundation to perform an extensive parameter fitting/fine-tuning and benchmarking framework for these models.

Therefore we present in the following all reported parameters and important data for comparing the models and their parameters in the future based on what the research community developed and benchmarked up to now.

COSMO-RS parameter tables ^{5 6}

As documented by Klamt et al. 1998 ⁵

Element-specific parameters

Element \(k\)	Cavity radius \(R_k\) \([Å]\)	Dispersion constant \(\gamma_k\) \([kcal/(mol*Å^2)]\)
H	1.30	-0.041
C	2.00	-0.037
N	1.83	-0.027
O	1.72	-0.042
Cl	2.05	-0.052

General COSMO-RS parameters

Parameter	Unit	Value
\(r_{av}\)	[Å]	0.5
\(a^´\)	[kcal/(mol Å^2)/e^2]	1288
\(f_{corr}\)	-	2.4
\(c_{hb}\)	[kcal/(mol Å^2/e^2)]	7400
\(\sigma_{hb}\)	[e/Å^2]	0.0082
\(a_{eff}\)	[Å^2]	7.1
\(\lambda\)	-	0.14
\(\omega\)	[kcal/mol]	-0.21
\(\nu\)	-	9.15

As documented by Gerlach et al. 2022 ⁶

Parameter	Unit	Turbomole	Orca-A	Orca-B
Basis set	-	TZVP	TZVPD	TZVPD
\(\sigma^{}\) applied	-	yes	no	yes
\(r_{av}\)	\([Å]\)	0.5 5	0.5 5	0.55
\(a_{eff}\)	\([Å^2]\)	6.25 7	5.943	6.115
\(\alpha_{mf}\)	\([kJ Å^2/(mol*e^2)]\)	5.950e6 7	5.953e6	7.584e6
\(f_{corr}\)	-	2.4 5		2.4
\(c_{hb}\)	\([kJ Å^2/(mol*e^2)]\)	3.670e7 7	3.671e7	3.093e7
\(\sigma_{hb}\)	\([e/Å^2]\)	8.5e-3 7	7.469e-3	7.876e-3
\(c_{hb^T}\)	-	1.5 8	1.5 8	1.5 8
\(A_{std}\)	\([Å^2]\)	79.53 9	51.28	41.89

openCOSMO-RS parameter tables ^{6 10 11}

As documented by Grigorash et al. 2025¹¹

Parameter	Unit	Value	References
\(a_{eff}\)	\([Å^2]\)	6.25 7
\(a_{eff}\)	\([Å^2]\)	6.226 []	ORCA
\(r_{av}\)	\([Å]\)	0.5 8
\(r_{av}^{cor}\)	\([Å]\)	1 5
\(\alpha\)	\([\frac{kJ*Å^2}{e^2*mol}]\)	5950.0e3 7	Turbomole
\(\alpha\)	\([\frac{kJ*Å^2}{e^2*mol}]\)	7.579075e6	ORCA
\(f_{corr}\)	-	2.4 5
\(\sigma_{hb}\)	\([\frac{e}{Å^2}]\)	8.5e-03 7	Turbomole
\(\sigma_{hb}\)	\([\frac{e}{Å^2}]\)	7.686e-03	ORCA
\(f_{e_I}\)	-	1 5
\(c_{hb}\)	\([\frac{kJ*Å^2}{e^2*mol}]\)	36700.0e3 7
\(c_{hb}\)	\([\frac{kJ*Å^2}{e^2*mol}]\)	2.7488747e7 7
\(c_{hb}^T\)	-	1.5 8
\(A_{std}\)	\([Å]\)	79.53 9	Turbomole
\(A_{std}\)	\([Å]\)	47.999	ORCA
\(z\)	-	10
\(a_{sge}\)	-	0.75	Turbomole

Parameterization of openCOSMO-RS 24a ¹⁰; gas and CPCM geometry optimizations at DFT/BP86/def2-TZVP level, gas and CPCM single point calculations at DFT/BP86/def2-TZVPD level with ORCA 6.0.

Parameter	Unit	Value
\(r_{av}^*\)	\([Å]\)	0.5
\(a_{eff}\)	\([Å^2]\)	5.925
\(\alpha_{mf}\)	\([kJ*Å/(mol*e^2)]\)	7281
\(f_{corr}^*\)	\([-]\)	2.4
\(c_{hb}\)	\([kJ Å^2/(mol*e^2)]\)	43327
\(c_{hb}^{T*}\)	\([-]\)	1.5
\(\sigma_{hb}\)	\([e/Å^2]\)	0.00961
\(A_{std}\)	\([Å^2]\)	41.624
\(\nu\)	\([kJ/mol]\)	-18.61
\(\omega_{ring}\)	\([kJ/mol]\)	1.100
\(\tau_1\)	\([kJ/(mol*Å^2)]\)	0.123
\(\tau_2\)	\([kJ/(mol*Å^2)]\)	0.096
\(\tau_3\)	\([kJ/(mol*Å^2)]\)	0.003
\(\tau_4\)	\([kJ/(mol*Å^2)]\)	0.015
\(\tau_5\)	\([kJ/(mol*Å^2)]\)	0.023
\(\tau_6\)	\([kJ/(mol*Å^2)]\)	0.143
\(\tau_7\)	\([kJ/(mol*Å^2)]\)	0.171
\(\tau_8\)	\([kJ/(mol*Å^2)]\)	0.018
\(\tau_9\)	\([kJ/(mol*Å^2)]\)	0.015
\(\tau_10\)	\([kJ/(mol*Å^2)]\)	0.146

Calculated parameters based on openCOSMO-RS implementation ⁶

Calculated parameters	Unit	Values
COSMOSPACE thresh		10e-06
COSMOSPACE maxit		1000
\(\sigma_{min} / \sigma^{\perp}_{min}\)		-0.15
\(\sigma_{max} / \sigma^{\perp}_{max}\)		0.15
\(\sigma_step / \sigma^{\perp}_{step}\)		0.001
HB donors		[1,100:150]
HB acceptors		[6,7,8,9,15,16,17,35,53]

Dispersion term support in openCOSMO-RS :

As documented by Grigorash et al. 2025

Combinatorial term	Unit	SG_6	SG_7	FH_6	FH_7	Elbro_6	Elbro_7
\(a_{eff}\)	\([Å^2]\)	6.115	6.115	5.034	5.034	2.745619	2.773896
\(r^*_{av}\)	\([Å^2]\)				0.5
\(\alpha_{mf}\)	\([J/(mol Å^2)/e^2]\)	7.548e+06	7.548e+06	7.592e+06	7.592e+06	1.210e+07	1.198e+07
\(f^*_{corr}\)					2.4
\(c_{hb}\)	\([J/(mol Å^2)/e^2]\)	3.093e+07	3.093e+07	3.094e+07	3.094e+07	3.093e+07 []	3.093e+07 []
\(c_{hb}^{T^*}\)					1.5
\(\sigma_{hb}\)	\([e/Å^2]\)	0.007876	0.007876	0.007276	0.007276	0.007876 []	0.007876 []
\(A_{std}\)	\([Å]\)	41.89	41.89
\(\tau_{C(sp^2)}^{vdW}\)	\(J^{0.5}/Å\)	11.193	9.425	14.123	10.577	17.675	18.221
\(\tau_{C(sp^2)^{vdW}}\)	\(J^{0.5}/Å\)		10.235		11.480		19.054
\(\tau_{H}^{vdW}\)	\(J^{0.5}/Å\)	10.041	9.021	12.581	10.300	15.636	17.506
\(\tau_{F}^{vdW}\)	\(J^{0.5}/Å\)	3.240	1.977	5.319	2.522	8.914	10.260
\(\tau_{Cl}^{vdW}\)	\(J^{0.5}/Å\)	11.865	10.647	14.575	11.735	18.289	19.669
\(\tau_{Br}^{vdW}\)	\(J^{0.5}/Å\)	17.602	16.414	20.796	17.545	25.965	27.218
\(\tau_{I}^{vdW}\)	\(J^{0.5}/Å\)	19.578	18.236	22.823	20.031	26.716	28.244

Cross-interaction parameter values and deviations for the respective openCOSMO-RS dispersion parametrizations

Parameters	SG_6_cross	SG_7_cross	FH_6_cross	FH_7_cross	Elbro_6_cross	Elbro_7_c
\(k_{H-C(sp^3)}\)	0.18159	0.01833	0.24871	0.00159	0.15927	-0.81875
\(k_{H-C(sp^2)}\)		-0.11976		0.15274		-0.34582
\(k_{H-F}\)	0.06897	0.66191	0.09735	0.74790	0.07372	-0.17848
\(k_{H-Cl}\)	0.22346	0.25723	0.16732	0.14724	0.14310	-0.04194
\(k_{H-Br}\)	-0.10650	-0.24431	0.16923	0.16498	0.10780	-0.25552
\(k_{H-I}\)	0.56261	0.29858	0.29177	0.22584	0.36641	0.11027
\(k_{C(sp^3)-C(sp^2)}\)		0.13688		0.01047		-0.07518
\(k_{C(sp^3)-F}\)	0.25838	-0.73798	0.24347	-0.94343	0.07947	-0.63680
\(k_{C(sp^2)-F}\)		-0.38252		0.01254		-0.14392
\(k_{C(sp^3)-Cl}\)	-0.09500	-0.19281	-0.01799	-0.10530	-0.03320	-0.33904
\(k_{C(sp^2)-Cl}\)		-0.19363		-0.03107		-0.11235
\(k_{C(sp^3)-Br}\)	-0.13948	0.07020	-0.08836	-0.20372	-0.09358	-0.04681
\(k_{C(sp^2)-Br}\)		-0.15348		0.00222		-0.04118
\(k_{C(sp^3)-I}\)	-0.47451	-0.25369	-0.25303	-0.21984	-0.30310	-0.45880
\(k_{C(sp^2)-I}\)		-0.42763		-0.37158		-0.28746
\(k_{F-Cl}\)	-0.02230	-0.14832	-0.03033	-0.10208	-0.01986	-0.02610
\(k_{F-Br}\)	-0.40404	-0.33715	-0.01564	-0.31896	-0.11410	0.02610
\(k_{F-I}\)	-0.12331	0.02270	-0.19483	0.05881	-0.12197	-0.11127
\(k_{Cl-Br}\)	0.11674	0.00928	-0.00186	0.06092	0.04101	-0.02610
\(k_{Cl-I}\)	-0.03422	-0.03802	-0.0192	-0.02616	-0.03001	-0.02558
\(k_{Br-I}\)	0.08028	-0.04237	0.01017	0.08432	0.05495	0.02172

Atomic polarizabilities ¹²

Element	Polarizability [a.u.]
H	4.5
C	11.3
F	3.74
Cl	14.6
Br	21
I	32.9

In a follow up publication Grigorash et al. introduced a new dispersion term based on atomic polarizabilities in a perfect conductor calculated as atomic polarizability tensors projected on atom segments.

The parameter set comparison documented by Grigorash et al. is shown in the following table:

Parameters	Unit	no dispersion	\(w_1\)	\(w_5\)	\(w_5\) + Eq. 20	\(w_5\) + Eq. 21
\(a_{eff}\)	\([Å^2]\)	4.574	4.943	4.955	4.90825	4.706
\(r_{av}*\)	\([Å]\)			0.5
\(\alpha_{mf}\)	\([kJ/(mol*Å^2*)/e^2]\)	7613	7294	7672	7876	7322
\(f_{corr}*\)	\([-]\)			2.4
\(c_{hb}\)	\([kJ/(mol*Å^2)/e^2]\)	53447	43265	43457	49318	43421
\(c_{hb}^{T*}\)	\([-]\)			1.5
\(\sigma_{hb}\)	\([e/Å^2]\)	0.009739	0.009492	0.00960	0.009953	0.009355
\(E_{corr}^F\)	\([J/mol]\)		349.18	350.37	346.82	340.73
\(m_{vdW}\)	\([-]\)		339.58	2.4097	29.567	27.853

Parameters	Unit	no dispersion	\(w_1\)	\(w_5\)	\(w_5\) + Eq. 20	\(w_5\) + Eq. 21
\(\tau_{1}\)	\([kJ/(mol*Å^2)]\)
\(\tau_{6}\)
\(\tau_{7}\)
\(\tau_{8}\)
\(\tau_{9}\)
\(\tau_{17}\)
\(\tau_{35}\)
\(\tau_{53}\)
\(\tau_{14}\)
\(\tau_{15}\)
\(\tau_{16}\)
\(\mu\)	\([kJ/mol]\)
\(\omega_{ring}\)	\([kJ/mol]\)

COSMO segment activity coefficient model: COSMO-SAC (2010)

The original COSMO-SAC model from Lin and Sandler ⁹ was published in 2002 and refined by Wang, Sandler and Chen in 2007 followed by further improvement by Hsieh, Sandler and Lin in 2010 ³. Conceptually and with respect to most model equations the COSMO-SAC is almost identical to the COSMO-RS model. We will highlight differences but most of the model performance deviations between COSMO-RS and COSMO-SAC are due to extensions and parameterizations with different experimental data sets.

In the following we will discuss the model equations for COSMO-SAC(2010) which has been further developed based on the COSMO-SAC(2007) model which has been refined with respect to: - Decoupled segments assumption and averaging methods - Definition of hydrogen bonding and the three classes of compounds - Regression of thermodynamic property data for pure compounds - Regression of thermodynamic vapour-liquid equilibrium data of binary and ternary mixtures

For more information on COSMO-SAC(2007), we recommend the respective paper .

The Gibbs energy is a foundational extensive thermodynamic state variable and the solvation Gibbs energy is made of summands which can be individually calculated via different models derived from quantum mechanical/chemical theory: $$ \Delta G_{sol} = \Delta G_{el} + \Delta G_{cav} + \Delta G_{disp} + \Delta G_{vib} + \Delta G_{lib} + \Delta G_{other} $$

To calculate the free solvation energy the COSMO-RS/-SAC models can be applied to quantify the charge density distribution of the cavity of a molecular structure. The COSMO-SAC(2007) model applied an empirical averaging method which is has been taken over by COSMO-SAC(2010):

\[\sigma_m = \frac{\sum_n \sigma_n^* (\frac{r_n^2 r_{eff}^2}{r_n^2 + r_{eff}^2}) exp[- f_{decay}(\frac{d_{mn}^2}{r_n^2+r_{eff}^2})]}{\sum_n \frac{r_n^2 r_{eff}^2}{r_n^2 + r_{eff}^2} + exp[-f_{decay}(\frac{d_{mn}^2}{r_n^2 + r_{eff}^2})]}\]

\(d_{mn}\) is the distance between segment \(m\) and \(n\). \(r_{eff}\) is the effective radius and in the COSMO-SAC(2010) publication ³ is mentioned that this should be the radius of the contact/standard surface area (\(r_{eff}=\sqrt{a_{eff}/\pi}\)).

An energy shift has to be applied because of the specific averaging method to correct for the dielectric change:

\[\Delta G_i*^{cc} = f_{pol}^{1/2}[E_{diel}(q) - E_{diel}(q*)]\]

where \(E_{diel}(q*) = 0.5 \sum_v \phi_v q_v*\) is the screening charge distribution and responsible for the dielectric polarization.

"In COSMO-RS and the COSMO-SAC(2002) model, the hydrogen-bonding effect in the segment exchange energy expression was included when the charge density of a segment exceeded an arbitrary threshold value \(\sigma_{hb}\) as follows" and Lin et al. have suggested to ... with the following hydrogen-bonding model equation by changing the COSMO-SAC(2002) hydrogen binding model:

\[ \begin{split} \Delta W(\sigma_m,\sigma_n) = &f_{pol} (\frac{0.3 a_{eff}^{3/2}}{2\epsilon_0})(\sigma_m + \sigma_n)^2 \\ &+ c_{hb}f_{hb}(T) min(0,min(\sigma_m,\sigma_n)+\sigma_{hb} max(0,max(\sigma_m,\sigma_n)-\sigma_{hb})) \end{split} \]

to

\[ \begin{split} \Delta W(\sigma_m^t, \sigma_n^s) = & f_{pol} (\frac{0.3 a_{eff}^{3/2}}{2 \epsilon} (\sigma_m^t + \sigma_n^s)^2) \\ & - c_{hb}(\sigma_m^t,\sigma_n^s)(\sigma_m^t - \sigma_n^s)^2 \end{split} \]

The superscripts \(t\) and \(s\) refer to a hydrogen bond (hb) or non-hydrogen bond (nhb) segment.

The \(c_{hb}(\sigma_m^t,\sigma_n^s)\) parameter is calculated as follows:

\[ c_{hb}(\sigma_m^t,\sigma_n^s) = \begin{cases} c_{hb} \text{ if } s=t=hb, \sigma_m^t * \sigma_n^s < 0 \\ 0 \end{cases} \]

The equation for calculating the activity coefficient with the COSMO-SAC(2007) model is also different from the COSMO-SAC(2002) model and has been changed from:

\[ ln \Gamma_{sol}(\sigma_m) = -ln( \sum_{\sigma_n} p(\sigma_n) \Gamma_{sol}(\sigma_n) exp(\frac{-\Delta W(\sigma_m,\sigma_n)}{RT})) \]

to

\[ ln(\Gamma_{sol}^t(\sigma_m^t)) = -ln(\sum_s^{nhb,hb} \sum_{\sigma_n} p^s(\sigma_n^s) \Gamma_{sol}^s(\sigma_n^s) exp(\frac{-\Delta W(\sigma_m^t,\sigma_n^s)}{RT})) \]

And then the solvation free energy can be calculated with the sigma profile data generated from the DFT-COSMO calculations for each molecule in a given mixture. The following equation from COSMO-SAC(2002):

\[ \frac{\Delta G*_{i/sol}^{res}}{RT} = n \sum_{\sigma_n} p(\sigma_m) ln(\Gamma_{sol}(\sigma_m)) \]

was modified to:

\[ \frac{\Delta G*_{i/sol}^{res}}{RT} = n \sum_{s}^{nhb,hb} p(\sigma_m) p(\sigma_m^s) ln(\Gamma_sol^s(\sigma_m^s)) \]

for COSMO-SAC(2007)

The COSMO-SAC(2010) introduced two modifications to the COSMO-SAC(2007) model:

• (I) The electrostatic interaction parameter cES is treated as a temperature-dependent parameter: $$ C_{ES} = A_{ES} + \frac{B_{ES}}{T^2} $$

• (II) The incorporation of hydrogen-bonding interactions for molecules where the hydrogen bonds with specific atom elements can be calculated correctly and thus differentiated: \(\(p_i^{hb}(\sigma) = p_i^{OH} + p_i^{OT}(\sigma)\)\)

with \(p_i^{OH}(\sigma)=(A_i^{OH}(\sigma)/A_i) x P^{HB}(\sigma)\) and \(p_i^{nhb}(\sigma) = (A^{OT}_i(\sigma)/A_i) x P^{HB}(\sigma)\)

With this second modification the surface type interactions are defined as:

\[ c_{hb}(\sigma_m^t, \sigma_n^s) = \begin{cases} c_{OH-OH} \text{ if } s=t=OH \text{ and } \sigma_m^t * \sigma_n^s < 0\\ c_{OT-OT} \text{ if } s=t=OH \text{ and } \sigma_m^t * \sigma_n^s < 0\\ c_{OH-OT} \text{ if } s=OH, t=OT, \text{ and } \sigma_m^t * \sigma_n^s <0\\ 0 \text{ otherwise } \end{cases} \]

Different parameters have been applied in the past, some even inconsistent and the parameter tables at the end of this COSMO-SAC section gives an overview of the parameter sets reported for COSMO-SAC(2007) and COSMO-SAC(2010).

COSMO-SAC (2007) and COSMO-SAC (2010) parameter tables

https://pubs.acs.org/doi/10.1021/ie404410v

Universal parameters

Parameter	Unit	COSMO-SAC (2007)	COSMO-SAC (2010) 3
\(a_{eff}\)	\(Å^2\)	7.25	7.25
\(f_{decay}\)	-	3.579	3.579
\(r_{avg}\)	\(Å\)	-	-
\(c_{es}\)	\(kcal/mol*Å^4/e^2\)	-	\(6326 + [1.5*10^8]/T^2\)
\(c_{hb}\)	\(kcal/mol*Å^4/e^2\)	3484	-
\(c_{OH-OH}\)	\(kcal/mol*Å^4/e^2\)	-	4014
\(c_{OH-OT}\)	\(kcal/mol*Å^4/e^2\)	-	3016
\(c_{OT-OT}\)	\(kcal/mol*Å^4/e^2\)	-	932
\(q_S\)	-	0.5	-
\(\sigma_0\)	\(e/Å^2\)	0.007	0.007
\(q_0\)	\(Å^2\)	79.53	79.53
\(r_0\)	\(Å^3\)	66.69	66.69

Atom bonding specific parameters and dispersion parameters (\(\epsilon/R^{\alpha}\) [\(K Å^3\)])

Atom type	\(R_i^{el}\) [Å]	COMSO-SAC (2007)	COSMO-SAC 2010 []
H	1.30	0
C (sp3)	2.00	36053.60
C (sp2)	2.00	30089.74
C (sp)	2.00	22133.70
N (sp3)	1.83	11859.65
N (sp2)	1.83	12693.51
N (sp)	1.83	2589.76
O(sp3,-H)	1.72	6759.19
O(sp3)	1.72	8420.93
O(sp2)	1.72	3289.98
O(sp2,-N)	1.72	10092.63
F	1.72	6741.91
P	2.12	100952.74 or 83439 []
S	2.16	61073.65 or 60949 []
Cl	2.05	35624.08
Br	2.16	60949.17 or 61074 []
I	2.32	83438.54 or 100953 []

---

1. Andreas Klamt, Gerard J. P. Krooshof, and Ross Taylor. Cosmospace: alternative to conventional activity-coefficient models. AIChE Journal, 48(10):2332–2349, 2002. doi:https://doi.org/10.1002/aic.690481023. ↩

2. A. J. Staverman. The entropy of high polymer solutions. generalization of formulae. Recueil des Travaux Chimiques des Pays-Bas, 69(2):163–174, 1950. doi:https://doi.org/10.1002/recl.19500690203. ↩

3. Chieh-Ming Hsieh, Stanley I. Sandler, and Shiang-Tai Lin. Improvements of cosmo-sac for vapor–liquid and liquid–liquid equilibrium predictions. Fluid Phase Equilibria, 297(1):90–97, 2010. doi:https://doi.org/10.1016/j.fluid.2010.06.011. ↩↩↩↩

4. Thomas Gerlach, Simon Müller, and Irina Smirnova. Development of a cosmo-rs based model for the calculation of phase equilibria in electrolyte systems. AIChE Journal, 64(1):272–285, 2018. doi:https://doi.org/10.1002/aic.15875. ↩

5. Andreas Klamt, Volker Jonas, Thorsten Bürger, and John C. W. Lohrenz. Refinement and parametrization of cosmo-rs. The Journal of Physical Chemistry A, 102(26):5074–5085, 1998. doi:10.1021/jp980017s. ↩↩↩↩↩↩↩↩↩

6. Thomas Gerlach, Simon Müller, Andrés González de Castilla, and Irina Smirnova. An open source cosmo-rs implementation and parameterization supporting the efficient implementation of multiple segment descriptors. Fluid Phase Equilibria, 560:113472, 2022. doi:https://doi.org/10.1016/j.fluid.2022.113472. ↩↩↩↩

7. Frank Eckert and Andreas Klamt. Fast solvent screening via quantum chemistry: cosmo-rs approach. AIChE Journal, 48(2):369–385, 2002. doi:https://doi.org/10.1002/aic.690480220. ↩↩↩↩↩↩↩↩↩

8. Andreas Klamt and Frank Eckert. Cosmo-rs: a novel and efficient method for the a priori prediction of thermophysical data of liquids. Fluid Phase Equilibria, 172(1):43–72, 2000. doi:https://doi.org/10.1016/S0378-3812(00)00357-5. ↩↩↩↩↩

9. Shiang-Tai Lin and Stanley I. Sandler. A priori phase equilibrium prediction from a segment contribution solvation model. Industrial & Engineering Chemistry Research, 41(5):899–913, 2002. doi:10.1021/ie001047w. ↩↩↩↩↩

10. Simon Müller, Thomas Nevolianis, Miquel Garcia-Ratés, Christoph Riplinger, Kai Leonhard, and Irina Smirnova. Predicting solvation free energies for neutral molecules in any solvent with opencosmo-rs. Fluid Phase Equilibria, 589:114250, 2025. doi:https://doi.org/10.1016/j.fluid.2024.114250. ↩↩

11. Daria Grigorash, Simon Müller, Patrice Paricaud, Erling H. Stenby, Irina Smirnova, and Wei Yan. A comprehensive approach to incorporating intermolecular dispersion into the opencosmo-rs model. part 1. halocarbons. Chemical Engineering Science, 309:121425, 2025. doi:https://doi.org/10.1016/j.ces.2025.121425. ↩↩

12. Peter Schwerdtfeger and Jeffrey K. Nagle and. 2018 table of static dipole polarizabilities of the neutral elements in the periodic table*. Molecular Physics, 117(9-12):1200–1225, 2019. doi:10.1080/00268976.2018.1535143. ↩

13. Shu Wang, Yuhua Song, and Chau-Chyun Chen. Extension of cosmo-sac solvation model for electrolytes. Industrial & Engineering Chemistry Research, 50(1):176–187, 2011. doi:10.1021/ie100689g. ↩

14. Thomas Ingram, Thomas Gerlach, Tanja Mehling, and Irina Smirnova. Extension of cosmo-rs for monoatomic electrolytes: modeling of liquid–liquid equilibria in presence of salts. Fluid Phase Equilibria, 314:29–37, 2012. doi:https://doi.org/10.1016/j.fluid.2011.09.021. ↩