feos-org · prehner · Apr 14, 2023 · Mar 22, 2023 · Mar 24, 2023 · Mar 25, 2023
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Changed
+- Changed the internal implementation of the association contribution to accomodate more general association schemes. [#150](https://github.com/feos-org/feos/pull/150)
+- To comply with the new association implementation, the default values of `na` and `nb` are now `0` rather than `1`. Parameter files have been adapted accordingly. [#150](https://github.com/feos-org/feos/pull/150)
 
 ## [0.4.2] - 2023-03-20
 - Python only: Release the changes introduced in `feos-core` 0.4.1 and `feos-dft` 0.4.1.

diff --git a/docs/theory/models/FeOs_Association.png b/docs/theory/models/FeOs_Association.png
diff --git a/docs/theory/models/association.md b/docs/theory/models/association.md
@@ -0,0 +1,83 @@
+# Association
+The Helmholtz contribution due to short range attractive interaction ("association") in SAFT models can be conveniently expressed as
+
+$$\frac{A^\mathrm{assoc}}{k_\mathrm{B}TV}=\sum_\alpha\rho_\alpha\left(\ln X_\alpha-\frac{X_\alpha}{2}+\frac{1}{2}\right)$$
+
+Here, $\alpha$ is the index of all distinguishable **association sites** in the system. The site density $\rho_\alpha$ is the density of the segment or molecule on which the association site is located times the multiplicity of the site. The fraction of non-bonded association sites $X_\alpha$ is calculated implicitly from
+
+$$\frac{1}{X_\alpha}=1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}$$ (eqn:x_assoc)
+
+where $\Delta^{\alpha\beta}$ is the association strength between sites $\alpha$ and $\beta$. The exact expression for the association strength varies between different SAFT versions. We implement the expressions used in PC-SAFT but a generalization to other models is straightforward.
+
+The number of degrees of freedom in the association strength matrix $\Delta$ is vast and for practical purposes it is useful to reduce the number of non-zero elements in $\Delta$. In $\text{FeO}_\text{s}$ we use 3 kinds of association sites: $A$, $B$ and $C$. Association sites of kind $A$ only associate with $B$ and vice versa, whereas association sites of kind $C$ only associate with other sites of kind $C$. By sorting the sites by their kind, the entire $\Delta$ matrix can be reduced to two smaller matrices: $\Delta^{AB}$ and $\Delta^{CC}$.
+
+```{image} FeOs_Association.png
+:alt: Association strength matrix
+:class: bg-primary
+:width: 300px
+:align: center
+```
+
+In practice, the $AB$ association can be linked to hydrogen bonding sites. The $CC$ association is less widely used but implemented to ensure that all the association schemes defined in [Huang and Radosz 1990](https://pubs.acs.org/doi/10.1021/ie00107a014) are covered.
+
+## Calculation of the fraction of non-bonded association sites
+
+The algorithm to solve the fraction of non-bonded sites for general association schemes follows [Michelsen 2006](https://pubs.acs.org/doi/full/10.1021/ie060029x). The problem is rewritten as an optimization problem with gradients
+
+$$g_\alpha=\frac{1}{X_\alpha}-1-\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}$$
+
+The analytic Hessian is
+
+$$H_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha^2}-\rho_\beta\Delta^{\alpha\beta}$$
+
+with the Kronecker delta $\delta_{\alpha\beta}$. However, [Michelsen 2006](https://pubs.acs.org/doi/full/10.1021/ie060029x) shows that it is beneficial to use
+
+$$\hat{H}_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha}\left(1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}\right)-\rho_\beta\Delta^{\alpha\beta}$$
+
+instead. $X_\alpha$ can then be solved robustly using a Newton algorithm with a check that ensures that the values of $X_\alpha$ remain positive. With the split in $AB$ and $CC$ association the two kinds different versions could be solved separately from each other. This is currently not implemented, because use cases are rare to nonexistent and the benefit is small.
+
+A drastic improvement in performance, however, can be achieved by solving eq. {eq}`eqn:x_assoc` analytically for simple cases. If there is only one $A$ and one $B$ site the corresponding fractions of non-bonded association sites can be calculated from
+
+$$X_A=\frac{2}{\sqrt{\left(1+\left(\rho_A-\rho_B\right)\Delta^{AB}\right)^2+4\rho_B\Delta^{AB}}+1+\left(\rho_B-\rho_A\right)\Delta^{AB}}$$
+$$X_B=\frac{2}{\sqrt{\left(1+\left(\rho_A-\rho_B\right)\Delta^{AB}\right)^2+4\rho_B\Delta^{AB}}+1+\left(\rho_A-\rho_B\right)\Delta^{AB}}$$
+
+The specific form of the expressions (with the square root in the denominator) is particularly robust for cases where $\Delta$ and/or $\rho$ are small or even 0.
+
+Analogously, for a single $C$ site, the expression becomes
+
+$$X_C=\frac{2}{\sqrt{1+4\rho_C\Delta^{CC}}+1}$$
+
+
+## PC-SAFT expression for the association strength
+In $\text{FeO}_\text{s}$ every association site $\alpha$ is parametrized with the dimensionless association volume $\kappa^\alpha$ and the association energy $\varepsilon^\alpha$. The association strength between sites $\alpha$ and $\beta$ is then calculated from
+
+$$\Delta^{\alpha\beta}=\left(\frac{1}{1-\zeta_3}+\frac{\frac{3}{2}d_{ij}\zeta_2}{\left(1-\zeta_3\right)^2}+\frac{\frac{1}{2}d_{ij}^2\zeta_2^2}{\left(1-\zeta_3\right)^3}\right)\sqrt{\sigma_i^3\kappa^\alpha\sigma_j^3\kappa^\beta}\left(e^{\frac{\varepsilon^\alpha+\varepsilon^\beta}{2k_\mathrm{B}T}}-1\right)$$
+
+with 
+
+$$d_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~d_k=\sigma_k\left(1-0.12e^{\frac{-3\varepsilon_k}{k_\mathrm{B}T}}\right)$$
+
+and
+
+$$\zeta_2=\frac{\pi}{6}\sum_k\rho_km_kd_k^2,~~~~\zeta_3=\frac{\pi}{6}\sum_k\rho_km_kd_k^3$$
+
+The indices $i$, $j$ and $k$ are used to index molecules or segments rather than sites. $i$ and $j$ in particular refer to the molecule or segment that contain site $\alpha$ or $\beta$ respectively.
+
+On their own the parameters $\kappa^\alpha$ and $\varepsilon^\beta$ have no physical meaning. For a pure system of self-associating molecules it is therefore natural to define $\kappa^A=\kappa^B\equiv\kappa^{A_iB_i}$ and $\varepsilon^A=\varepsilon^B\equiv\varepsilon^{A_iB_i}$ where $\kappa^{A_iB_i}$ and $\varepsilon^{A_iB_i}$ are now parameters describing the molecule rather than individual association sites. However, for systems that are not self-associating, the more generic parametrization is required.
+
+## Helmholtz energy functional
+The Helmholtz energy contribution proposed by [Yu and Wu 2002](https://aip.scitation.org/doi/abs/10.1063/1.1463435) is used to model association in inhomogeneous systems. It uses the same weighted densities that are used in [Fundamental Measure Theory](hard_spheres) (the White-Bear version). The Helmholtz energy density is then calculated from
+
+$$\beta f=\sum_\alpha\frac{n_0^\alpha\xi_i}{m_i}\left(\ln X_\alpha-\frac{X_\alpha}{2}+\frac{1}{2}\right)$$
+
+with the equations for the fraction of non-bonded association sites adapted to
+
+$$\frac{1}{X_\alpha}=1+\sum_\beta\frac{n_0^\beta\xi_j}{m_j}X_\beta\Delta^{\alpha\beta}$$
+
+The association strength, again using the PC-SAFT parametrization, is
+
+$$\Delta^{\alpha\beta}=\left(\frac{1}{1-n_3}+\frac{\frac{1}{4}d_{ij}n_2\xi}{\left(1-n_3\right)^2}+\frac{\frac{1}{72}d_{ij}^2n_2^2\xi}{\left(1-n_3\right)^3}\right)\sqrt{\sigma_i^3\kappa^\alpha\sigma_j^3\kappa^\beta}\left(e^{\frac{\varepsilon^\alpha+\varepsilon^\beta}{2k_\mathrm{B}T}}-1\right)$$
+
+The quantities $\xi$ and $\xi_i$ were introduced to account for the effect of inhomogeneity and are defined as
+
+$$\xi=1-\frac{\vec{n}_2\cdot\vec{n}_2}{n_2^2},~~~~\xi_i=1-\frac{\vec{n}_2^i\cdot\vec{n}_2^i}{{n_2^i}^2}$$
diff --git a/docs/theory/models/hard_spheres.md b/docs/theory/models/hard_spheres.md
@@ -0,0 +1,57 @@
+# Hard spheres
+
+$\text{FeO}_\text{s}$ provides an implementation of the Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) equation of state ([Boublík, 1970](https://doi.org/10.1063/1.1673824), [Mansoori et al., 1971](https://doi.org/10.1063/1.1675048)) for hard-sphere mixtures which is often used as reference contribution in SAFT equations of state. The implementation is generalized to allow the description of non-sperical or fused-sphere reference fluids.
+
+The reduced Helmholtz energy density is calculated according to
+
+$$\frac{\beta A}{V}=\frac{6}{\pi}\left(\frac{3\zeta_1\zeta_2}{1-\zeta_3}+\frac{\zeta_2^3}{\zeta_3\left(1-\zeta_3\right)^2}+\left(\frac{\zeta_2^3}{\zeta_3^2}-\zeta_0\right)\ln\left(1-\zeta_3\right)\right)$$
+
+with the packing fractions
+
+$$\zeta_k=\frac{\pi}{6}\sum_\alpha C_{k,\alpha}\rho_\alpha d_\alpha^k,~~~~~~~~k=0\ldots 3$$
+
+The geometry coefficients $C_{k,\alpha}$ and the segment diameters $d_\alpha$ depend on the context in which the model is used. The following table shows how the expression can be reused in various models. For details on the fused-sphere chain model, check the [repository](https://github.com/feos-org/feos-fused-chains) or the [publication](https://doi.org/10.1103/PhysRevE.105.034110).
+
+||Hard spheres|PC-SAFT|Fused-sphere chains|
+|-|:-:|:-:|:-:|
+|$d_\alpha$|$\sigma_\alpha$|$\sigma_\alpha\left(1-0.12e^{\frac{-3\varepsilon_\alpha}{k_\mathrm{B}T}}\right)$|$\sigma_\alpha$|
+|$C_{0,\alpha}$|$1$|$m_\alpha$|$1$|
+|$C_{1,\alpha}$|$1$|$m_\alpha$|$A_\alpha^*$|
+|$C_{1,\alpha}$|$1$|$m_\alpha$|$A_\alpha^*$|
+|$C_{1,\alpha}$|$1$|$m_\alpha$|$V_\alpha^*$|
+
+## Fundamental measure theory
+
+An model for inhomogeneous mixtues of hard spheres is provided by fundamental measure theory (FMT, [Rosenfeld, 1989](https://doi.org/10.1103/PhysRevLett.63.980)). Different variants have been proposed of which only those that are consistent with the BMCSL equation of state in the homogeneous limit are currently considered in $\text{FeO}_\text{s}$ (exluding, e.g., the original Rosenfeld and White-Bear II variants).
+
+The Helmholtz energy density is calculated according to
+
+$$\beta f=-n_0\ln\left(1-n_3\right)+\frac{n_{12}}{1-n_3}+\frac{1}{36\pi}n_2n_{22}f_3(n_3)$$
+
+The expressions for $n_{12}$ and $n_{22}$ depend on the FMT version.
+
+|version|$n_{12}$|$n_{22}$|references|
+|-|:-:|:-:|:-:|
+|WhiteBear|$n_1n_2-\vec n_1\cdot\vec n_2$|$n_2^2-3\vec n_2\cdot\vec n_2$|[Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313), [Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)|
+|KierlikRosinberg|$n_1n_2$|$n_2^2$|[Kierlik and Rosinberg, 1990](https://doi.org/10.1103/PhysRevA.42.3382)|
+|AntiSymWhiteBear|$n_1n_2-\vec n_1\cdot\vec n_2$|$n_2^2\left(1-\frac{\vec n_2\cdot\vec n_2}{n_2^2}\right)^3$|[Rosenfeld et al., 1997](https://doi.org/10.1103/PhysRevE.55.4245), [Kessler et al., 2021](https://doi.org/10.1016/j.micromeso.2021.111263)|
+
+For small $n_3$, the value of $f(n_3)$ numerically diverges. Therefore, it is approximated with a Taylor expansion.
+
+$$f_3=\begin{cases}\frac{n_3+\left(1-n_3\right)^2\ln\left(1-n_3\right)}{n_3^2\left(1-n_3\right)^2}&\text{if }n_3>10^{-5}\\
+\frac{3}{2}+\frac{8}{3}n_3+\frac{15}{4}n_3^2+\frac{24}{5}n_3^3+\frac{35}{6}n_3^4&\text{else}\end{cases}$$
+
+The weighted densities $n_k(\mathbf{r})$ are calculated by convolving the density profiles $\rho_\alpha(\mathbf{r})$ with weight functions $\omega_k^\alpha(\mathbf{r})$
+
+$$n_k(\mathbf{r})=\sum_\alpha n_k^\alpha(\mathbf{r})=\sum_\alpha\int\rho_\alpha(\mathbf{r}')\omega_k^\alpha(\mathbf{r}-\mathbf{r}')\mathrm{d}\mathbf{r}'$$
+
+which differ between the different FMT versions.
+
+||WhiteBear/AntiSymWhiteBear|KierlikRosinberg|
+|-|:-:|:-:|
+|$\omega_0^\alpha(\mathbf{r})$|$\frac{C_{0,\alpha}}{\pi\sigma_\alpha^2}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{0,\alpha}\left(-\frac{1}{8\pi}\,\delta''\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)+\frac{1}{2\pi\|\mathbf{r}\|}\,\delta'\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)\right)$|
+|$\omega_1^\alpha(\mathbf{r})$|$\frac{C_{1,\alpha}}{2\pi\sigma_\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$\frac{C_{1,\alpha}}{8\pi}\,\delta'\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|
+|$\omega_2^\alpha(\mathbf{r})$|$C_{2,\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{2,\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|
+|$\omega_3^\alpha(\mathbf{r})$|$C_{3,\alpha}\,\Theta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{3,\alpha}\,\Theta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|
+|$\vec\omega_1^\alpha(\mathbf{r})$|$C_{3,\alpha}\frac{\mathbf{r}}{2\pi\sigma_\alpha\|\mathbf{r}\|}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|-|
+|$\vec\omega_2^\alpha(\mathbf{r})$|$C_{3,\alpha}\frac{\mathbf{r}}{\|\mathbf{r}\|}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|-|
diff --git a/docs/theory/models/index.md b/docs/theory/models/index.md
@@ -6,4 +6,7 @@ It is currently still under construction. You can help by [contributing](https:/
 ```{eval-rst}
 .. toctree::
    :maxdepth: 1
+
+   hard_spheres
+   association
 ```
diff --git a/parameters/pcsaft/gross2002.json b/parameters/pcsaft/gross2002.json
@@ -13,7 +13,9 @@
             "sigma": 3.23,
             "epsilon_k": 188.9,
             "kappa_ab": 0.035176,
-            "epsilon_k_ab": 2899.5
+            "epsilon_k_ab": 2899.5,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 32.042
     },
@@ -31,7 +33,9 @@
             "sigma": 3.1771,
             "epsilon_k": 198.24,
             "kappa_ab": 0.032384,
-            "epsilon_k_ab": 2653.4
+            "epsilon_k_ab": 2653.4,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 46.069
     },
@@ -49,7 +53,9 @@
             "sigma": 3.2522,
             "epsilon_k": 233.4,
             "kappa_ab": 0.015268,
-            "epsilon_k_ab": 2276.8
+            "epsilon_k_ab": 2276.8,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 60.096
     },
@@ -67,7 +73,9 @@
             "sigma": 3.6139,
             "epsilon_k": 259.59,
             "kappa_ab": 0.006692,
-            "epsilon_k_ab": 2544.6
+            "epsilon_k_ab": 2544.6,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 74.123
     },
@@ -85,7 +93,9 @@
             "sigma": 3.4508,
             "epsilon_k": 247.28,
             "kappa_ab": 0.010319,
-            "epsilon_k_ab": 2252.1
+            "epsilon_k_ab": 2252.1,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 88.15
     },
@@ -103,7 +113,9 @@
             "sigma": 3.6735,
             "epsilon_k": 262.32,
             "kappa_ab": 0.005747,
-            "epsilon_k_ab": 2538.9
+            "epsilon_k_ab": 2538.9,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 102.177
     },
@@ -121,7 +133,9 @@
             "sigma": 3.545,
             "epsilon_k": 253.46,
             "kappa_ab": 0.001155,
-            "epsilon_k_ab": 2878.5
+            "epsilon_k_ab": 2878.5,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 116.203
     },
@@ -139,7 +153,9 @@
             "sigma": 3.7145,
             "epsilon_k": 262.74,
             "kappa_ab": 0.002197,
-            "epsilon_k_ab": 2754.8
+            "epsilon_k_ab": 2754.8,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 130.23
     },
@@ -157,7 +173,9 @@
             "sigma": 3.7292,
             "epsilon_k": 263.64,
             "kappa_ab": 0.001427,
-            "epsilon_k_ab": 2941.9
+            "epsilon_k_ab": 2941.9,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 144.257
     },
@@ -175,7 +193,9 @@
             "sigma": 3.2085,
             "epsilon_k": 208.42,
             "kappa_ab": 0.024675,
-            "epsilon_k_ab": 2253.9
+            "epsilon_k_ab": 2253.9,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 60.096
     },
@@ -193,7 +213,9 @@
             "sigma": 3.9053,
             "epsilon_k": 266.01,
             "kappa_ab": 0.001863,
-            "epsilon_k_ab": 2618.8
+            "epsilon_k_ab": 2618.8,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 88.15
     },
@@ -211,7 +233,9 @@
             "sigma": 3.0007,
             "epsilon_k": 366.51,
             "kappa_ab": 0.034868,
-            "epsilon_k_ab": 2500.7
+            "epsilon_k_ab": 2500.7,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 18.015
     },
@@ -229,7 +253,9 @@
             "sigma": 2.8906,
             "epsilon_k": 214.94,
             "kappa_ab": 0.095103,
-            "epsilon_k_ab": 684.3
+            "epsilon_k_ab": 684.3,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 31.06
     },
@@ -247,7 +273,9 @@
             "sigma": 3.1343,
             "epsilon_k": 221.53,
             "kappa_ab": 0.017275,
-            "epsilon_k_ab": 854.7
+            "epsilon_k_ab": 854.7,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 45.09
     },
@@ -265,7 +293,9 @@
             "sigma": 3.5347,
             "epsilon_k": 250.52,
             "kappa_ab": 0.022674,
-            "epsilon_k_ab": 1028.1
+            "epsilon_k_ab": 1028.1,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 59.11
     },
@@ -283,7 +313,9 @@
             "sigma": 3.4777,
             "epsilon_k": 231.8,
             "kappa_ab": 0.02134,
-            "epsilon_k_ab": 932.2
+            "epsilon_k_ab": 932.2,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 59.11
     },
@@ -301,7 +333,9 @@
             "sigma": 3.7021,
             "epsilon_k": 335.47,
             "kappa_ab": 0.074883,
-            "epsilon_k_ab": 1351.6
+            "epsilon_k_ab": 1351.6,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 93.13
     },
@@ -319,7 +353,9 @@
             "sigma": 3.8582,
             "epsilon_k": 211.59,
             "kappa_ab": 0.07555,
-            "epsilon_k_ab": 3044.4
+            "epsilon_k_ab": 3044.4,
+            "na": 1.0,
+            "nb": 1.0
         },
         "molarweight": 60.053
     }

diff --git a/parameters/pcsaft/loetgeringlin2015_homo.json b/parameters/pcsaft/loetgeringlin2015_homo.json
@@ -277,6 +277,8 @@
       "epsilon_k": 488.66,
       "epsilon_k_ab": 2517.0,
       "kappa_ab": 0.006825,
+      "na": 1.0,
+      "nb": 1.0,
       "viscosity": [
         -15.7583e-3,
         -2.5654e-1,
@@ -294,6 +296,8 @@
       "epsilon_k": 467.59,
       "epsilon_k_ab": 1064.6,
       "kappa_ab": 0.026662,
+      "na": 1.0,
+      "nb": 1.0,
       "viscosity": [
         -4.4048e-3,
         -0.6089e-1,