feos_pcsaft/dft/
pure_saft_functional.rs

1use super::association::N0_CUTOFF;
2use super::polar::{pair_integral_ij, triplet_integral_ijk};
3use crate::eos::association::{assoc_site_frac_a, assoc_site_frac_ab};
4use crate::eos::dispersion::{A0, A1, A2, B0, B1, B2};
5use crate::eos::polar::{AD, AQ, BD, BQ, CD, CQ, PI_SQ_43};
6use crate::parameters::PcSaftParameters;
7use feos_core::{EosError, EosResult};
8use feos_dft::fundamental_measure_theory::FMTVersion;
9use feos_dft::{
10    FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape,
11};
12use ndarray::*;
13use num_dual::*;
14use std::f64::consts::{FRAC_PI_6, PI};
15use std::fmt;
16use std::rc::Rc;
17
18const PI36M1: f64 = 1.0 / (36.0 * PI);
19const N3_CUTOFF: f64 = 1e-5;
20
21#[derive(Clone)]
22pub struct PureFMTAssocFunctional {
23    parameters: Rc<PcSaftParameters>,
24    version: FMTVersion,
25}
26
27impl PureFMTAssocFunctional {
28    pub fn new(parameters: Rc<PcSaftParameters>, version: FMTVersion) -> Self {
29        Self {
30            parameters,
31            version,
32        }
33    }
34}
35
36impl<N: DualNum<f64> + ScalarOperand> FunctionalContributionDual<N> for PureFMTAssocFunctional {
37    fn weight_functions(&self, temperature: N) -> WeightFunctionInfo<N> {
38        let r = self.parameters.hs_diameter(temperature) * 0.5;
39        WeightFunctionInfo::new(arr1(&[0]), false).extend(
40            vec![
41                WeightFunctionShape::Delta,
42                WeightFunctionShape::Theta,
43                WeightFunctionShape::DeltaVec,
44            ]
45            .into_iter()
46            .map(|s| WeightFunction {
47                prefactor: self.parameters.m.mapv(|m| m.into()),
48                kernel_radius: r.clone(),
49                shape: s,
50            })
51            .collect(),
52            false,
53        )
54    }
55
56    fn calculate_helmholtz_energy_density(
57        &self,
58        temperature: N,
59        weighted_densities: ArrayView2<N>,
60    ) -> EosResult<Array1<N>> {
61        let p = &self.parameters;
62
63        // weighted densities
64        let n2 = weighted_densities.index_axis(Axis(0), 0);
65        let n3 = weighted_densities.index_axis(Axis(0), 1);
66        let n2v = weighted_densities.slice_axis(Axis(0), Slice::new(2, None, 1));
67
68        // temperature dependent segment radius
69        let r = self.parameters.hs_diameter(temperature)[0] * 0.5;
70
71        // auxiliary variables
72        if n3.iter().any(|n3| n3.re() > 1.0) {
73            return Err(EosError::IterationFailed(String::from(
74                "PureFMTAssocFunctional",
75            )));
76        }
77        let ln31 = n3.mapv(|n3| (-n3).ln_1p());
78        let n3rec = n3.mapv(|n3| n3.recip());
79        let n3m1 = n3.mapv(|n3| -n3 + 1.0);
80        let n3m1rec = n3m1.mapv(|n3m1| n3m1.recip());
81        let n1 = n2.mapv(|n2| n2 / (r * 4.0 * PI));
82        let n0 = n2.mapv(|n2| n2 / (r * r * 4.0 * PI));
83        let n1v = n2v.mapv(|n2v| n2v / (r * 4.0 * PI));
84
85        let (n1n2, n2n2) = match self.version {
86            FMTVersion::WhiteBear => (
87                &n1 * &n2 - (&n1v * &n2v).sum_axis(Axis(0)),
88                &n2 * &n2 - (&n2v * &n2v).sum_axis(Axis(0)) * 3.0,
89            ),
90            FMTVersion::AntiSymWhiteBear => {
91                let mut xi2 = (&n2v * &n2v).sum_axis(Axis(0)) / n2.map(|n| n.powi(2));
92                xi2.iter_mut().for_each(|x| {
93                    if x.re() > 1.0 {
94                        *x = N::one()
95                    }
96                });
97                (
98                    &n1 * &n2 - (&n1v * &n2v).sum_axis(Axis(0)),
99                    &n2 * &n2 * xi2.mapv(|x| (-x + 1.0).powi(3)),
100                )
101            }
102            _ => unreachable!(),
103        };
104
105        // The f3 term contains a 0/0, therefore a taylor expansion is used for small values of n3
106        let mut f3 = (&n3m1 * &n3m1 * &ln31 + n3) * &n3rec * n3rec * &n3m1rec * &n3m1rec;
107        f3.iter_mut().zip(n3).for_each(|(f3, &n3)| {
108            if n3.re() < N3_CUTOFF {
109                *f3 = (((n3 * 35.0 / 6.0 + 4.8) * n3 + 3.75) * n3 + 8.0 / 3.0) * n3 + 1.5;
110            }
111        });
112        let mut phi = -(&n0 * &ln31) + n1n2 * &n3m1rec + n2n2 * n2 * PI36M1 * f3;
113
114        // association
115        if p.nassoc == 1 {
116            let mut xi = -(&n2v * &n2v).sum_axis(Axis(0)) / (&n2 * &n2) + 1.0;
117            xi.iter_mut().zip(&n2).for_each(|(xi, &n2)| {
118                if n2.re() < N0_CUTOFF * 4.0 * PI * p.m[0] * r.re().powi(2) {
119                    *xi = N::one();
120                }
121            });
122
123            let k = &n2 * &n3m1rec * r;
124            let deltarho = (((&k / 18.0 + 0.5) * &k * &xi + 1.0) * n3m1rec)
125                * ((temperature.recip() * p.epsilon_k_aibj[(0, 0)]).exp_m1()
126                    * (p.sigma[0].powi(3) * p.kappa_aibj[(0, 0)]))
127                * (&n0 / p.m[0] * &xi);
128
129            let f = |x: N| x.ln() - x * 0.5 + 0.5;
130            phi = phi
131                + if p.nb[0] > 0.0 {
132                    let xa = deltarho.mapv(|d| assoc_site_frac_ab(d, p.na[0], p.nb[0]));
133                    let xb = (xa.clone() - 1.0) * p.na[0] / p.nb[0] + 1.0;
134                    (n0 / p.m[0] * xi) * (xa.mapv(f) * p.na[0] + xb.mapv(f) * p.nb[0])
135                } else {
136                    let xa = deltarho.mapv(|d| assoc_site_frac_a(d, p.na[0]));
137                    n0 / p.m[0] * xi * (xa.mapv(f) * p.na[0])
138                };
139        }
140
141        Ok(phi)
142    }
143}
144
145impl fmt::Display for PureFMTAssocFunctional {
146    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
147        write!(f, "Pure FMT+association")
148    }
149}
150
151#[derive(Clone)]
152pub struct PureChainFunctional {
153    parameters: Rc<PcSaftParameters>,
154}
155
156impl PureChainFunctional {
157    pub fn new(parameters: Rc<PcSaftParameters>) -> Self {
158        Self { parameters }
159    }
160}
161
162impl<N: DualNum<f64> + ScalarOperand> FunctionalContributionDual<N> for PureChainFunctional {
163    fn weight_functions(&self, temperature: N) -> WeightFunctionInfo<N> {
164        let d = self.parameters.hs_diameter(temperature);
165        WeightFunctionInfo::new(arr1(&[0]), true)
166            .add(
167                WeightFunction::new_scaled(d.clone(), WeightFunctionShape::Delta),
168                false,
169            )
170            .add(
171                WeightFunction {
172                    prefactor: (&self.parameters.m / 8.0).mapv(|x| x.into()),
173                    kernel_radius: d,
174                    shape: WeightFunctionShape::Theta,
175                },
176                false,
177            )
178    }
179
180    fn calculate_helmholtz_energy_density(
181        &self,
182        _: N,
183        weighted_densities: ArrayView2<N>,
184    ) -> EosResult<Array1<N>> {
185        let rho = weighted_densities.index_axis(Axis(0), 0);
186        // negative lambdas lead to nan, therefore the absolute value is used
187        let lambda = weighted_densities
188            .index_axis(Axis(0), 1)
189            .map(|&l| if l.re() < 0.0 { -l } else { l } + N::from(f64::EPSILON));
190        let eta = weighted_densities.index_axis(Axis(0), 2);
191
192        let y = eta.mapv(|eta| (eta * 0.5 - 1.0) / (eta - 1.0).powi(3));
193        Ok(-(y * lambda).mapv(|x| (x.ln() - 1.0) * (self.parameters.m[0] - 1.0)) * rho)
194    }
195}
196
197impl fmt::Display for PureChainFunctional {
198    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
199        write!(f, "Pure chain")
200    }
201}
202
203#[derive(Clone)]
204pub struct PureAttFunctional {
205    parameters: Rc<PcSaftParameters>,
206}
207
208impl PureAttFunctional {
209    pub fn new(parameters: Rc<PcSaftParameters>) -> Self {
210        Self { parameters }
211    }
212}
213
214impl<N: DualNum<f64> + ScalarOperand> FunctionalContributionDual<N> for PureAttFunctional {
215    fn weight_functions(&self, temperature: N) -> WeightFunctionInfo<N> {
216        let d = self.parameters.hs_diameter(temperature);
217        const PSI: f64 = 1.3862; // Homosegmented DFT (Sauer2017)
218        WeightFunctionInfo::new(arr1(&[0]), false).add(
219            WeightFunction::new_scaled(d * PSI, WeightFunctionShape::Theta),
220            false,
221        )
222    }
223
224    fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo<N> {
225        let d = self.parameters.hs_diameter(temperature);
226        const PSI: f64 = 1.3286; // pDGT (Rehner2018)
227        WeightFunctionInfo::new(arr1(&[0]), false).add(
228            WeightFunction::new_scaled(d * PSI, WeightFunctionShape::Theta),
229            false,
230        )
231    }
232
233    fn calculate_helmholtz_energy_density(
234        &self,
235        temperature: N,
236        weighted_densities: ArrayView2<N>,
237    ) -> EosResult<Array1<N>> {
238        let p = &self.parameters;
239        let rho = weighted_densities.index_axis(Axis(0), 0);
240
241        // temperature dependent segment radius
242        let d = p.hs_diameter(temperature)[0];
243
244        let eta = rho.mapv(|rho| rho * FRAC_PI_6 * p.m[0] * d.powi(3));
245        let m1 = (p.m[0] - 1.0) / p.m[0];
246        let m2 = m1 * (p.m[0] - 2.0) / p.m[0];
247        let e = temperature.recip() * p.epsilon_k[0];
248        let s3 = p.sigma[0].powi(3);
249
250        // I1, I2 and C1
251        let mut i1: Array1<N> = Array::zeros(eta.raw_dim());
252        let mut i2: Array1<N> = Array::zeros(eta.raw_dim());
253        for i in 0..=6 {
254            i1 = i1 + eta.mapv(|eta| eta.powi(i as i32) * (A0[i] + m1 * A1[i] + m2 * A2[i]));
255            i2 = i2 + eta.mapv(|eta| eta.powi(i as i32) * (B0[i] + m1 * B1[i] + m2 * B2[i]));
256        }
257        let c1 = eta.mapv(|eta| {
258            ((eta * 8.0 - eta.powi(2) * 2.0) / (eta - 1.0).powi(4) * p.m[0]
259                + (eta * 20.0 - eta.powi(2) * 27.0 + eta.powi(3) * 12.0 - eta.powi(4) * 2.0)
260                    / ((eta - 1.0) * (eta - 2.0)).powi(2)
261                    * (1.0 - p.m[0])
262                + 1.0)
263                .recip()
264        });
265        let mut phi = rho.mapv(|rho| -(rho * p.m[0]).powi(2) * e * s3 * PI)
266            * (i1 * 2.0 + c1 * i2.mapv(|i2| i2 * p.m[0] * e));
267
268        // dipoles
269        if p.ndipole > 0 {
270            let mu2_term = e * s3 * p.mu2[0];
271            let m = p.m[0].min(2.0);
272            let m1 = (m - 1.0) / m;
273            let m2 = m1 * (m - 2.0) / m;
274
275            let phi2 = -(&rho * &rho)
276                * pair_integral_ij(m1, m2, &eta, &AD, &BD, e)
277                * (mu2_term * mu2_term / s3 * PI);
278            let phi3 = -(&rho * &rho * rho)
279                * triplet_integral_ijk(m1, m2, &eta, &CD)
280                * (mu2_term * mu2_term * mu2_term / s3 * PI_SQ_43);
281
282            let mut phi_d = &phi2 * &phi2 / (&phi2 - &phi3);
283            phi_d.iter_mut().zip(phi2.iter()).for_each(|(p, &p2)| {
284                if p.re().is_nan() {
285                    *p = p2;
286                }
287            });
288            phi += &phi_d;
289        }
290
291        // quadrupoles
292        if p.nquadpole > 0 {
293            let q2_term = e * p.sigma[0].powi(5) * p.q2[0];
294            let m = p.m[0].min(2.0);
295            let m1 = (m - 1.0) / m;
296            let m2 = m1 * (m - 2.0) / m;
297
298            let phi2 = -(&rho * &rho)
299                * pair_integral_ij(m1, m2, &eta, &AQ, &BQ, e)
300                * (q2_term * q2_term / p.sigma[0].powi(7) * PI * 0.5625);
301            let phi3 = (&rho * &rho * rho)
302                * triplet_integral_ijk(m1, m2, &eta, &CQ)
303                * (q2_term * q2_term * q2_term / s3.powi(3) * PI * PI * 0.5625);
304
305            let mut phi_q = &phi2 * &phi2 / (&phi2 - &phi3);
306            phi_q.iter_mut().zip(phi2.iter()).for_each(|(p, &p2)| {
307                if p.re().is_nan() {
308                    *p = p2;
309                }
310            });
311            phi += &phi_q;
312        }
313
314        Ok(phi)
315    }
316}
317
318impl fmt::Display for PureAttFunctional {
319    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
320        write!(f, "Pure attractive")
321    }
322}
feos_pcsaft/dft/pure_saft_functional.rs

feos_pcsaft/dft/
pure_saft_functional.rs
