feos 0.9.5

FeOs - A framework for equations of state and classical density functional theory.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

use super::polar::helmholtz_energy_density_polar;
use crate::hard_sphere::HardSphereProperties;
use crate::pcsaft::eos::dispersion::{A0, A1, A2, B0, B1, B2};
use crate::pcsaft::parameters::PcSaftPars;
use feos_core::FeosError;
use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape};
use nalgebra::DVector;
use ndarray::*;
use num_dual::DualNum;
use std::f64::consts::{FRAC_PI_6, PI};

/// psi Parameter for DFT (Sauer2017)
const PSI_DFT: f64 = 1.3862;
/// psi Parameter for pDGT (Rehner2018)
const PSI_PDGT: f64 = 1.3286;

#[derive(Clone)]
pub(crate) struct AttractiveFunctional<'a> {
    parameters: &'a PcSaftPars,
}

impl<'a> AttractiveFunctional<'a> {
    pub(crate) fn new(parameters: &'a PcSaftPars) -> Self {
        Self { parameters }
    }
}

fn att_weight_functions<N: DualNum<f64> + Copy>(
    p: &PcSaftPars,
    psi: f64,
    temperature: N,
) -> WeightFunctionInfo<N> {
    let d = p.hs_diameter(temperature);
    let psi = N::from(psi);
    WeightFunctionInfo::new(DVector::from_fn(d.len(), |i, _| i), false).add(
        WeightFunction::new_scaled(d * psi, WeightFunctionShape::Theta),
        false,
    )
}

impl<'a> FunctionalContribution for AttractiveFunctional<'a> {
    fn name(&self) -> &'static str {
        "Attractive functional"
    }

    fn weight_functions<N: DualNum<f64> + Copy>(&self, temperature: N) -> WeightFunctionInfo<N> {
        att_weight_functions(self.parameters, PSI_DFT, temperature)
    }

    fn weight_functions_pdgt<N: DualNum<f64> + Copy>(
        &self,
        temperature: N,
    ) -> WeightFunctionInfo<N> {
        att_weight_functions(self.parameters, PSI_PDGT, temperature)
    }

    fn helmholtz_energy_density<N: DualNum<f64> + Copy>(
        &self,
        temperature: N,
        density: ArrayView2<N>,
    ) -> Result<Array1<N>, FeosError> {
        // auxiliary variables
        let p = &self.parameters;
        let n = p.m.len();

        // temperature dependent segment radius
        let d = p.hs_diameter(temperature);

        // packing fraction
        let d3m = d.zip_map(&p.m, |d, m| d.powi(3) * (m * FRAC_PI_6));
        let eta = density.outer_iter().zip(&d3m).fold(
            Array::zeros(density.raw_dim().remove_axis(Axis(0))),
            |acc: Array1<N>, (rho, &d3m)| acc + &rho * d3m,
        );

        // mean segment number
        let mut rhog = Array::zeros(eta.raw_dim());
        let mut m_bar = Array::zeros(eta.raw_dim());
        for (rhoi, &mi) in density.axis_iter(Axis(0)).zip(p.m.iter()) {
            m_bar += &(&rhoi * mi);
            rhog += &rhoi;
        }
        m_bar.iter_mut().zip(rhog.iter()).for_each(|(m, &r)| {
            if r.re() > f64::EPSILON {
                *m /= r
            } else {
                *m = N::one()
            }
        });

        // mixture densities, crosswise interactions of all segments on all chains
        let mut rho1mix: Array1<N> = Array::zeros(eta.raw_dim());
        let mut rho2mix: Array1<N> = Array::zeros(eta.raw_dim());
        for i in 0..n {
            for j in 0..n {
                let eps_ij_t = temperature.recip() * p.epsilon_k_ij[(i, j)];
                let sigma_ij_3 = p.sigma_ij[(i, j)].powi(3);
                rho1mix = rho1mix
                    + (&density.index_axis(Axis(0), i) * &density.index_axis(Axis(0), j))
                        .mapv(|x| x * (eps_ij_t * sigma_ij_3 * p.m[i] * p.m[j]));
                rho2mix = rho2mix
                    + (&density.index_axis(Axis(0), i) * &density.index_axis(Axis(0), j))
                        .mapv(|x| x * (eps_ij_t * eps_ij_t * sigma_ij_3 * p.m[i] * p.m[j]));
            }
        }

        // I1, I2 and C1
        let mut i1: Array1<N> = Array::zeros(eta.raw_dim());
        let mut i2: Array1<N> = Array::zeros(eta.raw_dim());
        let mut eta_i: Array1<N> = Array::ones(eta.raw_dim());
        let m1 = (m_bar.clone() - 1.0) / &m_bar;
        let m2 = (m_bar.clone() - 2.0) / &m_bar * &m1;
        for i in 0..=6 {
            i1 = i1 + (&m2 * A2[i] + &m1 * A1[i] + A0[i]) * &eta_i;
            i2 = i2 + (&m2 * B2[i] + &m1 * B1[i] + B0[i]) * &eta_i;
            eta_i = &eta_i * &eta;
        }
        let c1 = Zip::from(&eta).and(&m_bar).map_collect(|&eta, &m| {
            (m * (eta * 8.0 - eta.powi(2) * 2.0) / (eta - 1.0).powi(4)
                + (eta * (eta * (eta * (eta * 2.0 - 12.0) + 27.0) - 20.0))
                    / ((eta - 1.0) * (eta - 2.0)).powi(2)
                    * (m - 1.0)
                + 1.0)
                .recip()
        });

        // Helmholtz energy density
        let phi_polar = helmholtz_energy_density_polar(p, temperature, density)?;
        Ok((-rho1mix * i1 * 2.0 - rho2mix * m_bar * c1 * i2) * PI + phi_polar)
    }
}