COSMO-SAC Activity Coefficient Model in MATLAB

This repository contains a MATLAB implementation of the COSMO-SAC (COSMO Segment Activity Coefficient) model. The code calculates the thermodynamic activity coefficients of components in chemical mixtures based on their molecular surface charge distributions, known as $\sigma$-profiles.

---

Overview

The model predicts the non-ideal behavior of liquid mixtures by analyzing the interaction energies between molecular surface segments. It accounts for:

• Electrostatic Misfit Energy: Interactions between segments with different charge densities.
• Hydrogen Bonding: Specific attractions between highly polar surface regions based on a cutoff value.
• Combinatorial Effects: Differences in molecular size and shape using the Staverman-Guggenheim term.

---

Mathematical Logic

The calculation follows a multi-step thermodynamic workflow derived from the source code:

1. Sigma Profile Calculation ($\sigma$-profile)

For each molecule, the code calculates the surface area distribution as a function of charge density $\sigma$. The local segment charge is averaged over a circular area to account for neighbor effects using the following logic:

$$\sigma_{new, i} = \frac{\sum_{j} \sigma_j \frac{RAD_j^2 REFF^2}{RAD_j^2 + REFF^2} \exp \left( \frac{-Term1^2}{RAD_j^2 + REFF^2} \right)}{\sum_{j} \frac{RAD_j^2 REFF^2 (RAD_j^2 + REFF^2)}{RAD_j^2 + REFF^2} \exp \left( \frac{-Term1^2}{RAD_j^2 + REFF^2} \right)}$$

Where:

• Term1 is the Euclidean distance between surface segments $i$ and $j$.
• REFF is the effective averaging radius ($0.817642$ Å).

2. Interaction Energy ($\Delta w$)

The exchange energy between two segments $i$ and $j$ is calculated as:

$$\Delta w_{ij} = \frac{\alpha'}{2}(\sigma_i + \sigma_j)^2 + c_{hb} \max(0, \sigma_{acc} - \sigma_{hb}) \min(0, \sigma_{don} + \sigma_{hb})$$

• Misfit term: Determined by ALPHAPRIME and the sum of charge densities.
• HB term: Active when charge densities exceed the threshold SIGMAHB ($0.0084$ e/Ų).

3. Segment Activity Coefficients ($\Gamma_s$)

The activity coefficient of a segment is solved iteratively until the change is less than $1 \times 10^{-6}$:

$$\ln \Gamma_s(\sigma_i) = -\ln \left[ \sum_{j} P(\sigma_j) \Gamma_s(\sigma_j) \exp \left( \frac{-\Delta w_{ij}}{RT} \right) \right]$$

The code computes this for both pure species (SEGGAMMAPR) and the mixture (SEGGAMMA).

4. Total Activity Coefficient ($\gamma$)

The final activity coefficient for component $n$ is the exponential sum of the residual and combinatorial parts:

$$\gamma_n = \exp(\ln \gamma_n^{res} + \ln \gamma_n^{sg})$$

• Residual: Calculated by summing the difference between mixture and pure segment activity coefficients weighted by the sigma profile.
• Combinatorial (SG): Uses COORD (coordination number), RNORM (normalized volume), and QNORM (normalized area) parameters.

---

File Structure

File	Description
eqCOSMO.m	Main script for Two Phase Equilibrium; defines temperature and components.
Binary.m	Primary function calculating mole fractions, pure/mixture segment gammas, and final $\gamma$.
SimgaProfileCalculator.m	Manages the averaging of surface charges and calls the sorting logic.
Library.m	Reads geometry, volume, and charge data from inputQM.xlsx.
paraCOSMO.m	Provides universal constants ($R$, $c_{hb}$, $\sigma_{hb}$, $\alpha$, etc.).
SortSimgaProfile.m	Bins segment data into a discrete density distribution from -0.025 to 0.025.
ConvertAU2A.m	Converts coordinates from Atomic Units to Angstroms.

---

Usage Instructions

Prerequisites

• MATLAB installed.
• An Excel file named inputQM.xlsx in the working directory.
• The Excel file must have sheets named after your components (e.g., '67-63-0-2d' or 'Z1') containing the required QM output data.

Execution

1. Open eqCOSMO.m.
2. Set your desired SYSTEMP (Temperature in K).
3. Update ListCOMP with the sheet names from your Excel file.
4. Run the script.

Outputs

• sProfiles.xlsx: Stores the calculated $\sigma$-densities and profiles for each component.
• MixGamma.xlsx: Stores the final matrix containing mole fractions ($x$), activity coefficients ($\gamma$), and $\ln \gamma$.