Skip to content

Commit 1dc165a

Browse files
g-bauerprehner
g-bauer
and
prehner
authored
Implement SAFT-VR Mie equation of state (#237)
* PC SAFT done. Checking benches * Changes according to review, added SAFTVRQMie * adjust dft * Python fixes * Add uvtheory to Cargo feature * Started impl of SAFT VR Mie * Working dispersion and chain. Association needs evaluation. * Cleanup, removal of old files * Improved diameter calculation * Finished impl. Added test cases. * Silence warning of unused method * Fixed bug from rebasing * Fixed benchmark * Updated Readme * Updated changelog * Added parameters of Lafitte et al 2013, readme and bib file * Fixed integration test * Added example notebook with validation of implementation * Added docs * Changed visibility in core and dft. Fixed bug in benches. * Added section for association term in documentation * Added more examples to notebook * removed unused method from hard sphere module, added tests to CI, moved vr mie benches, removed fmt from Cargo.toml --------- Co-authored-by: Philipp Rehner <prehner@ethz.ch>
1 parent a19fbf4 commit 1dc165a

File tree

29 files changed

+4118
-51
lines changed

29 files changed

+4118
-51
lines changed

‎.github/workflows/test.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ jobs:
2828
strategy:
2929
fail-fast: false
3030
matrix:
31-
model: [pcsaft, gc_pcsaft, pets, uvtheory, saftvrqmie]
31+
model: [pcsaft, gc_pcsaft, pets, uvtheory, saftvrqmie, saftvrmie]
3232

3333
steps:
3434
- uses: actions/checkout@v4

‎CHANGELOG.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
55
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
66

77
## [Unreleased]
8+
### Added
9+
- Added SAFT-VR Mie equation of state.[#237](https://github.com/feos-org/feos/pull/237)
810

911
### Changed
1012
- 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‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,9 +66,10 @@ gc_pcsaft = ["association"]
6666
uvtheory = ["lazy_static"]
6767
pets = []
6868
saftvrqmie = []
69+
saftvrmie = []
6970
rayon = ["dep:rayon", "ndarray/rayon", "feos-core/rayon", "feos-dft?/rayon"]
7071
python = ["pyo3", "numpy", "quantity/python", "feos-core/python", "feos-dft?/python", "rayon"]
71-
all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie"]
72+
all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets", "saftvrqmie", "saftvrmie"]
7273

7374
[[bench]]
7475
name = "state_properties"
@@ -82,6 +83,10 @@ harness = false
8283
name = "dual_numbers"
8384
harness = false
8485

86+
[[bench]]
87+
name = "dual_numbers_saftvrmie"
88+
harness = false
89+
8590
[[bench]]
8691
name = "contributions"
8792
harness = false

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ The following models are currently published as part of the `FeOs` framework
4040
|`pets`|perturbed truncated and shifted Lennard-Jones mixtures|✓|✓|
4141
|`uvtheory`|equation of state for Mie fluids and mixtures|✓||
4242
|`saftvrqmie`|equation of state for quantum fluids and mixtures|✓|✓|
43+
|`saftvrmie`|statistical associating fluid theory for variable range interactions of Mie form|✓||
4344

4445
The list is being expanded continuously. Currently under development are implementations of ePC-SAFT and a Helmholtz energy functional for the UV theory.
4546

‎benches/dual_numbers_saftvrmie.rs‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
//! Benchmarks for the evaluation of the Helmholtz energy function
2+
//! for a given `StateHD` for different types of dual numbers.
3+
//! These should give an idea about the expected slow-down depending
4+
//! on the dual number type used without the overhead of the `State`
5+
//! creation.
6+
use criterion::{criterion_group, criterion_main, Criterion};
7+
use feos::hard_sphere::HardSphereProperties;
8+
use feos::saftvrmie::{test_utils::test_parameters, SaftVRMie, SaftVRMieParameters};
9+
use feos_core::si::*;
10+
use feos_core::{Derivative, Residual, State, StateHD};
11+
use ndarray::{Array, ScalarOperand};
12+
use num_dual::DualNum;
13+
use std::sync::Arc;
14+
15+
/// Helper function to create a state for given parameters.
16+
/// - temperature is 80% of critical temperature,
17+
/// - volume is critical volume,
18+
/// - molefracs (or moles) for equimolar mixture.
19+
fn state_saftvrmie(parameters: &Arc<SaftVRMieParameters>) -> State<SaftVRMie> {
20+
let n = parameters.pure_records.len();
21+
let eos = Arc::new(SaftVRMie::new(parameters.clone()));
22+
let moles = Array::from_elem(n, 1.0 / n as f64) * 10.0 * MOL;
23+
let cp = State::critical_point(&eos, Some(&moles), None, Default::default()).unwrap();
24+
let temperature = 0.8 * cp.temperature;
25+
State::new_nvt(&eos, temperature, cp.volume, &moles).unwrap()
26+
}
27+
28+
/// Residual Helmholtz energy given an equation of state and a StateHD.
29+
fn a_res<D: DualNum<f64> + Copy + ScalarOperand, E: Residual>(inp: (&Arc<E>, &StateHD<D>)) -> D {
30+
inp.0.residual_helmholtz_energy(inp.1)
31+
}
32+
33+
fn d_hs<D: DualNum<f64> + Copy>(inp: (&SaftVRMieParameters, D)) -> D {
34+
inp.0.hs_diameter(inp.1)[0]
35+
}
36+
37+
/// Benchmark for evaluation of the Helmholtz energy for different dual number types.
38+
fn bench_dual_numbers<E: Residual>(
39+
c: &mut Criterion,
40+
group_name: &str,
41+
state: State<E>,
42+
parameters: &SaftVRMieParameters,
43+
) {
44+
let mut group = c.benchmark_group(group_name);
45+
group.bench_function("d_f64", |b| {
46+
b.iter(|| d_hs((&parameters, state.derive0().temperature)))
47+
});
48+
group.bench_function("d_dual", |b| {
49+
b.iter(|| d_hs((&parameters, state.derive1(Derivative::DT).temperature)))
50+
});
51+
group.bench_function("d_dual2", |b| {
52+
b.iter(|| d_hs((&parameters, state.derive2(Derivative::DT).temperature)))
53+
});
54+
group.bench_function("d_hyperdual", |b| {
55+
b.iter(|| {
56+
d_hs((
57+
&parameters,
58+
state
59+
.derive2_mixed(Derivative::DT, Derivative::DT)
60+
.temperature,
61+
))
62+
})
63+
});
64+
group.bench_function("d_dual3", |b| {
65+
b.iter(|| d_hs((&parameters, state.derive3(Derivative::DT).temperature)))
66+
});
67+
68+
group.bench_function("a_f64", |b| {
69+
b.iter(|| a_res((&state.eos, &state.derive0())))
70+
});
71+
group.bench_function("a_dual", |b| {
72+
b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV))))
73+
});
74+
group.bench_function("a_dual2", |b| {
75+
b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV))))
76+
});
77+
group.bench_function("a_hyperdual", |b| {
78+
b.iter(|| {
79+
a_res((
80+
&state.eos,
81+
&state.derive2_mixed(Derivative::DV, Derivative::DV),
82+
))
83+
})
84+
});
85+
group.bench_function("a_dual3", |b| {
86+
b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV))))
87+
});
88+
}
89+
90+
/// Benchmark for the SAFT VR Mie equation of state
91+
fn saftvrmie(c: &mut Criterion) {
92+
let parameters = Arc::new(test_parameters().remove("ethane").unwrap());
93+
bench_dual_numbers(
94+
c,
95+
"dual_numbers_saftvrmie_ethane",
96+
state_saftvrmie(&parameters),
97+
&parameters,
98+
);
99+
}
100+
101+
criterion_group!(bench, saftvrmie);
102+
criterion_main!(bench);

‎docs/api/eos.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ If you want to adjust parameters of a model to experimental data you can use cla
2121
EquationOfState.python_residual
2222
EquationOfState.python_ideal_gas
2323
EquationOfState.uvtheory
24+
EquationOfState.saftvrmie
2425
EquationOfState.saftvrqmie
2526
```
2627

@@ -54,7 +55,7 @@ If you want to adjust parameters of a model to experimental data you can use cla
5455

5556
## The `estimator` module
5657

57-
### Import
58+
### Import
5859

5960
```python
6061
from feos.eos.estimator import Estimator, DataSet, Loss, Phase

‎docs/api/index.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,8 @@ These modules contain the objects to e.g. read parameters from files or build pa
2424
peng_robinson
2525
pets
2626
uvtheory
27+
saftvrmie
2728
saftvrqmie
2829
joback
2930
dippr
30-
```
31+
```

‎docs/api/saftvrmie.md‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
# `feos.saftvrmie`
2+
3+
Utilities to build `SaftVRMieParameters`.
4+
5+
## Example
6+
7+
```python
8+
from feos.saftvrmie import SaftVRMieParameters
9+
10+
path = 'parameters/saftvrmie/lafitte2013.json'
11+
parameters = SaftVRMieParameters.from_json(['ethane', 'methane'], path)
12+
```
13+
14+
## Data types
15+
16+
```{eval-rst}
17+
.. currentmodule:: feos.saftvrmie
18+
19+
.. autosummary::
20+
:toctree: generated/
21+
22+
Identifier
23+
ChemicalRecord
24+
PureRecord
25+
BinaryRecord
26+
SaftVRMieRecord
27+
SaftVRMieParameters
28+
```

‎docs/theory/models/association.md‎

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
# Association
2+
23
The Helmholtz contribution due to short range attractive interaction ("association") in SAFT models can be conveniently expressed as
34

45
$$\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
4748

4849
$$X_C=\frac{2}{\sqrt{1+4\rho_C\Delta^{CC}}+1}$$
4950

50-
5151
## PC-SAFT expression for the association strength
52+
5253
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
5354

5455
$$\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)$$
5556

56-
with
57+
with
5758

5859
$$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)$$
5960

@@ -65,7 +66,31 @@ The indices $i$, $j$ and $k$ are used to index molecules or segments rather than
6566

6667
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.
6768

69+
## SAFT-VR Mie expression for the association strength
70+
71+
We provide an implementation of the association strength as published by [Lafitte et al. (2013)](https://doi.org/10.1063/1.4819786).
72+
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.
73+
74+
The association strength between sites $\alpha$ and $\beta$ is then calculated from
75+
76+
$$\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)$$
77+
78+
with
79+
80+
$$\tilde{d}_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}$$
81+
82+
and
83+
84+
$$\zeta_2=\frac{\pi}{6}\sum_k\rho_km_kd_k^2,~~~~\zeta_3=\frac{\pi}{6}\sum_k\rho_km_kd_k^3$$
85+
86+
where
87+
88+
$$d_i = \int_{l}^{\sigma_i} \left[1 - e^{-\beta u_i^\text{Mie}(r)}\right]\mathrm{d}r.$$
89+
90+
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.
91+
6892
## Helmholtz energy functional
93+
6994
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
7095

7196
$$\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
80105

81106
The quantities $\xi$ and $\xi_i$ were introduced to account for the effect of inhomogeneity and are defined as
82107

83-
$$\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}$$
108+
$$\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}$$

0 commit comments

Comments
 (0)