Implement SAFT-VR Mie equation of state

.github/workflows/test.yml

@@ -28,7 +28,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        model: [pcsaft, gc_pcsaft, pets, uvtheory, saftvrqmie]
+        model: [pcsaft, gc_pcsaft, pets, uvtheory, saftvrqmie, saftvrmie]
 
     steps:
       - uses: actions/checkout@v4

CHANGELOG.md

@@ -5,6 +5,8 @@ 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]
+### Added
+- Added SAFT-VR Mie equation of state.[#237](https://github.com/feos-org/feos/pull/237)
 
 ### Changed
 - Updated model implementations to account for the removal of trait objects for Helmholtz energy contributions and the de Broglie in `feos-core`. [#226](https://github.com/feos-org/feos/pull/226)

Cargo.toml

@@ -66,9 +66,10 @@ gc_pcsaft = ["association"]
 uvtheory = ["lazy_static"]
 pets = []
 saftvrqmie = []
+saftvrmie = []
 rayon = ["dep:rayon", "ndarray/rayon", "feos-core/rayon", "feos-dft?/rayon"]
 python = ["pyo3", "numpy", "quantity/python", "feos-core/python", "feos-dft?/python", "rayon"]
-all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie"]
+all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie", "saftvrmie"]
 
 [[bench]]
 name = "state_properties"
@@ -82,6 +83,10 @@ harness = false
 name = "dual_numbers"
 harness = false
 
+[[bench]]
+name = "dual_numbers_saftvrmie"
+harness = false
+
 [[bench]]
 name = "contributions"
 harness = false

README.md

@@ -40,6 +40,7 @@ The following models are currently published as part of the `FeOs` framework
 |`pets`|perturbed truncated and shifted Lennard-Jones mixtures|✓|✓|
 |`uvtheory`|equation of state for Mie fluids and mixtures|✓||
 |`saftvrqmie`|equation of state for quantum fluids and mixtures|✓|✓|
+|`saftvrmie`|statistical associating fluid theory for variable range interactions of Mie form|✓||
 
 The list is being expanded continuously. Currently under development are implementations of ePC-SAFT and a Helmholtz energy functional for the UV theory.
 

benches/dual_numbers_saftvrmie.rs

@@ -0,0 +1,102 @@
+//! Benchmarks for the evaluation of the Helmholtz energy function
+//! for a given `StateHD` for different types of dual numbers.
+//! These should give an idea about the expected slow-down depending
+//! on the dual number type used without the overhead of the `State`
+//! creation.
+use criterion::{criterion_group, criterion_main, Criterion};
+use feos::hard_sphere::HardSphereProperties;
+use feos::saftvrmie::{test_utils::test_parameters, SaftVRMie, SaftVRMieParameters};
+use feos_core::si::*;
+use feos_core::{Derivative, Residual, State, StateHD};
+use ndarray::{Array, ScalarOperand};
+use num_dual::DualNum;
+use std::sync::Arc;
+
+/// Helper function to create a state for given parameters.
+/// - temperature is 80% of critical temperature,
+/// - volume is critical volume,
+/// - molefracs (or moles) for equimolar mixture.
+fn state_saftvrmie(parameters: &Arc<SaftVRMieParameters>) -> State<SaftVRMie> {
+    let n = parameters.pure_records.len();
+    let eos = Arc::new(SaftVRMie::new(parameters.clone()));
+    let moles = Array::from_elem(n, 1.0 / n as f64) * 10.0 * MOL;
+    let cp = State::critical_point(&eos, Some(&moles), None, Default::default()).unwrap();
+    let temperature = 0.8 * cp.temperature;
+    State::new_nvt(&eos, temperature, cp.volume, &moles).unwrap()
+}
+
+/// Residual Helmholtz energy given an equation of state and a StateHD.
+fn a_res<D: DualNum<f64> + Copy + ScalarOperand, E: Residual>(inp: (&Arc<E>, &StateHD<D>)) -> D {
+    inp.0.residual_helmholtz_energy(inp.1)
+}
+
+fn d_hs<D: DualNum<f64> + Copy>(inp: (&SaftVRMieParameters, D)) -> D {
+    inp.0.hs_diameter(inp.1)[0]
+}
+
+/// Benchmark for evaluation of the Helmholtz energy for different dual number types.
+fn bench_dual_numbers<E: Residual>(
+    c: &mut Criterion,
+    group_name: &str,
+    state: State<E>,
+    parameters: &SaftVRMieParameters,
+) {
+    let mut group = c.benchmark_group(group_name);
+    group.bench_function("d_f64", |b| {
+        b.iter(|| d_hs((&parameters, state.derive0().temperature)))
+    });
+    group.bench_function("d_dual", |b| {
+        b.iter(|| d_hs((&parameters, state.derive1(Derivative::DT).temperature)))
+    });
+    group.bench_function("d_dual2", |b| {
+        b.iter(|| d_hs((&parameters, state.derive2(Derivative::DT).temperature)))
+    });
+    group.bench_function("d_hyperdual", |b| {
+        b.iter(|| {
+            d_hs((
+                &parameters,
+                state
+                    .derive2_mixed(Derivative::DT, Derivative::DT)
+                    .temperature,
+            ))
+        })
+    });
+    group.bench_function("d_dual3", |b| {
+        b.iter(|| d_hs((&parameters, state.derive3(Derivative::DT).temperature)))
+    });
+
+    group.bench_function("a_f64", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive0())))
+    });
+    group.bench_function("a_dual", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV))))
+    });
+    group.bench_function("a_dual2", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV))))
+    });
+    group.bench_function("a_hyperdual", |b| {
+        b.iter(|| {
+            a_res((
+                &state.eos,
+                &state.derive2_mixed(Derivative::DV, Derivative::DV),
+            ))
+        })
+    });
+    group.bench_function("a_dual3", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV))))
+    });
+}
+
+/// Benchmark for the SAFT VR Mie equation of state
+fn saftvrmie(c: &mut Criterion) {
+    let parameters = Arc::new(test_parameters().remove("ethane").unwrap());
+    bench_dual_numbers(
+        c,
+        "dual_numbers_saftvrmie_ethane",
+        state_saftvrmie(&parameters),
+        &parameters,
+    );
+}
+
+criterion_group!(bench, saftvrmie);
+criterion_main!(bench);

docs/api/eos.md

@@ -21,6 +21,7 @@ If you want to adjust parameters of a model to experimental data you can use cla
     EquationOfState.python_residual
     EquationOfState.python_ideal_gas
     EquationOfState.uvtheory
+    EquationOfState.saftvrmie
     EquationOfState.saftvrqmie
 ```
 
@@ -54,7 +55,7 @@ If you want to adjust parameters of a model to experimental data you can use cla
 
 ## The `estimator` module
 
-### Import 
+### Import
 
 ```python
 from feos.eos.estimator import Estimator, DataSet, Loss, Phase

docs/api/index.md

@@ -24,7 +24,8 @@ These modules contain the objects to e.g. read parameters from files or build pa
    peng_robinson
    pets
    uvtheory
+   saftvrmie
    saftvrqmie
    joback
    dippr
-```
+```

docs/api/saftvrmie.md

@@ -0,0 +1,28 @@
+# `feos.saftvrmie`
+
+Utilities to build `SaftVRMieParameters`.
+
+## Example
+
+```python
+from feos.saftvrmie import SaftVRMieParameters
+
+path = 'parameters/saftvrmie/lafitte2013.json'
+parameters = SaftVRMieParameters.from_json(['ethane', 'methane'], path)
+```
+
+## Data types
+
+```{eval-rst}
+.. currentmodule:: feos.saftvrmie
+
+.. autosummary::
+    :toctree: generated/
+
+    Identifier
+    ChemicalRecord
+    PureRecord
+    BinaryRecord
+    SaftVRMieRecord
+    SaftVRMieParameters
+```

docs/theory/models/association.md

@@ -1,4 +1,5 @@
 # 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)$$
@@ -47,13 +48,13 @@ 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 
+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)$$
 
@@ -65,7 +66,31 @@ The indices $i$, $j$ and $k$ are used to index molecules or segments rather than
 
 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.
 
+## SAFT-VR Mie expression for the association strength
+
+We provide an implementation of the association strength as published by [Lafitte et al. (2013)](https://doi.org/10.1063/1.4819786).
+Every association site $\alpha$ is parametrized with two distances $r_c^\alpha$ and $r_d^\alpha$, and the association energy $\varepsilon^\alpha$. Whereas $r_c^\alpha$ is a parameter that is adjusted for each substance, $r_d^\alpha$ is kept constant as $r_d^\alpha = 0.4 \sigma$. Note that the parameter $r_c^\alpha$ has to be provided as dimensionless quantity in the input, i.e. divided by the segment's $\sigma$ value. 
+
+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}\tilde{d}_{ij}\zeta_2}{\left(1-\zeta_3\right)^2}+\frac{\frac{1}{2}\tilde{d}_{ij}^2\zeta_2^2}{\left(1-\zeta_3\right)^3}\right)K_{ij}^{\alpha\beta}\sigma_{ij}^3\left(e^{\frac{\sqrt{\varepsilon^\alpha\varepsilon^\beta}}{k_\mathrm{B}T}}-1\right)$$
+
+with
+
+$$\tilde{d}_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}$$
+
+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$$
+
+where 
+
+$$d_i = \int_{l}^{\sigma_i} \left[1 - e^{-\beta u_i^\text{Mie}(r)}\right]\mathrm{d}r.$$
+
+The integral is solved using a Gauss-Legendre quadrature of fifth order where the lower limit, $l$, is determined using the method of [Aasen et al.](https://doi.org/10.1063/1.5111364) The dimensionless bonding volume $K^{\alpha\beta}_{ij}$ (see [A43](https://doi.org/10.1063/1.4819786)) utilizes arithmetic combining rules for $r_c^{\alpha\beta}$, $r_d^{\alpha\beta}$ and $d_{ij}$ of unlike sites and segments, respectively.
+
 ## 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)$$
@@ -80,4 +105,4 @@ $$\Delta^{\alpha\beta}=\left(\frac{1}{1-n_3}+\frac{\frac{1}{4}d_{ij}n_2\xi}{\lef
 
 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}$$
+$$\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}$$