feos-pets 0.1.0

Implementation of PeTS equation of state and corresponding Helmholtz energy functional.
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

use crate::eos::dispersion::{A, B};
use crate::parameters::PetsParameters;
use feos_core::EosError;
use feos_dft::{
    FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape,
};
use ndarray::*;
use num_dual::DualNum;
use std::f64::consts::{FRAC_PI_3, PI};
use std::fmt;
use std::rc::Rc;

/// psi Parameter for DFT (Heier2018)
const PSI_DFT: f64 = 1.21;
/// psi Parameter for pDGT (not adjusted, yet)
const PSI_PDGT: f64 = 1.21;

#[derive(Clone)]
pub struct AttractiveFunctional {
    parameters: Rc<PetsParameters>,
}

impl AttractiveFunctional {
    pub fn new(parameters: Rc<PetsParameters>) -> Self {
        Self {
            parameters: parameters.clone(),
        }
    }
}

fn att_weight_functions<N: DualNum<f64> + ScalarOperand>(
    p: &PetsParameters,
    psi: f64,
    temperature: N,
) -> WeightFunctionInfo<N> {
    let d = p.hs_diameter(temperature);
    WeightFunctionInfo::new(Array1::from_shape_fn(d.len(), |i| i), false).add(
        WeightFunction::new_scaled(d * psi, WeightFunctionShape::Theta),
        false,
    )
}

impl<N: DualNum<f64> + ScalarOperand> FunctionalContributionDual<N> for AttractiveFunctional {
    fn weight_functions(&self, temperature: N) -> WeightFunctionInfo<N> {
        att_weight_functions(&self.parameters, PSI_DFT, temperature)
    }

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

    fn calculate_helmholtz_energy_density(
        &self,
        temperature: N,
        density: ArrayView2<N>,
    ) -> Result<Array1<N>, EosError> {
        // auxiliary variables
        let p = &self.parameters;
        let n = p.sigma.len();

        // temperature dependent segment radius
        let r = p.hs_diameter(temperature) * 0.5;

        // packing fraction
        let eta = density
            .outer_iter()
            .zip((&r * &r * &r * 4.0 * FRAC_PI_3).into_iter())
            .fold(
                Array::zeros(density.raw_dim().remove_axis(Axis(0))),
                |acc: Array1<N>, (rho, r3)| acc + &rho * r3,
            );

        // 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));
                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));
            }
        }

        // 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());
        for i in 0..=6 {
            i1 = i1 + &eta_i * A[i];
            i2 = i2 + &eta_i * B[i];
            eta_i = &eta_i * &eta;
        }
        let c1 =
            eta.mapv(|eta| ((eta * 8.0 - eta.powi(2) * 2.0) / (eta - 1.0).powi(4) + 1.0).recip());

        // Helmholtz energy density
        Ok((-rho1mix * i1 * 2.0 - rho2mix * c1 * i2) * PI)
    }
}

impl fmt::Display for AttractiveFunctional {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "Attractive functional")
    }
}