feos_dft/
functional_contribution.rs

1use crate::weight_functions::WeightFunctionInfo;
2use feos_core::{EosResult, StateHD};
3use ndarray::prelude::*;
4use ndarray::{RemoveAxis, ScalarOperand};
5use num_dual::*;
6use num_traits::Zero;
7use std::fmt::Display;
8
9/// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers.
10pub trait FunctionalContribution: Display + Sync + Send {
11    /// Return the weight functions required in this contribution.
12    fn weight_functions<N: DualNum<f64> + Copy + ScalarOperand>(
13        &self,
14        temperature: N,
15    ) -> WeightFunctionInfo<N>;
16
17    /// Overwrite this if the weight functions in pDGT are different than for DFT.
18    fn weight_functions_pdgt<N: DualNum<f64> + Copy + ScalarOperand>(
19        &self,
20        temperature: N,
21    ) -> WeightFunctionInfo<N> {
22        self.weight_functions(temperature)
23    }
24
25    /// Return the Helmholtz energy density for the given temperature and weighted densities.
26    fn helmholtz_energy_density<N: DualNum<f64> + Copy + ScalarOperand>(
27        &self,
28        temperature: N,
29        weighted_densities: ArrayView2<N>,
30    ) -> EosResult<Array1<N>>;
31
32    fn helmholtz_energy<N: DualNum<f64> + Copy + ScalarOperand>(&self, state: &StateHD<N>) -> N {
33        // calculate weight functions
34        let weight_functions = self.weight_functions(state.temperature);
35
36        // calculate segment density
37        let density = weight_functions
38            .component_index
39            .mapv(|c| state.partial_density[c]);
40
41        // calculate weighted density and Helmholtz energy
42        let weight_constants = weight_functions.weight_constants(Zero::zero(), 0);
43        let weighted_densities = weight_constants.dot(&density).insert_axis(Axis(1));
44        self.helmholtz_energy_density(state.temperature, weighted_densities.view())
45            .unwrap()[0]
46            * state.volume
47    }
48
49    fn first_partial_derivatives<N: DualNum<f64> + Copy>(
50        &self,
51        temperature: N,
52        weighted_densities: Array2<N>,
53        mut helmholtz_energy_density: ArrayViewMut1<N>,
54        mut first_partial_derivative: ArrayViewMut2<N>,
55    ) -> EosResult<()> {
56        let mut wd = weighted_densities.mapv(Dual::from_re);
57        let t = Dual::from_re(temperature);
58        let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
59
60        for i in 0..wd.shape()[0] {
61            wd.index_axis_mut(Axis(0), i)
62                .map_inplace(|x| x.eps = N::one());
63            phi = self.helmholtz_energy_density(t, wd.view())?;
64            first_partial_derivative
65                .index_axis_mut(Axis(0), i)
66                .assign(&phi.mapv(|p| p.eps));
67            wd.index_axis_mut(Axis(0), i)
68                .map_inplace(|x| x.eps = N::zero());
69        }
70        helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
71        Ok(())
72    }
73
74    fn second_partial_derivatives(
75        &self,
76        temperature: f64,
77        weighted_densities: ArrayView2<f64>,
78        mut helmholtz_energy_density: ArrayViewMut1<f64>,
79        mut first_partial_derivative: ArrayViewMut2<f64>,
80        mut second_partial_derivative: ArrayViewMut3<f64>,
81    ) -> EosResult<()> {
82        let mut wd = weighted_densities.mapv(HyperDual64::from);
83        let t = HyperDual64::from(temperature);
84        let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
85
86        for i in 0..wd.shape()[0] {
87            wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps1 = 1.0);
88            for j in 0..=i {
89                wd.index_axis_mut(Axis(0), j).map_inplace(|x| x.eps2 = 1.0);
90                phi = self.helmholtz_energy_density(t, wd.view())?;
91                let p = phi.mapv(|p| p.eps1eps2);
92                second_partial_derivative
93                    .index_axis_mut(Axis(0), i)
94                    .index_axis_mut(Axis(0), j)
95                    .assign(&p);
96                if i != j {
97                    second_partial_derivative
98                        .index_axis_mut(Axis(0), j)
99                        .index_axis_mut(Axis(0), i)
100                        .assign(&p);
101                }
102                wd.index_axis_mut(Axis(0), j).map_inplace(|x| x.eps2 = 0.0);
103            }
104            first_partial_derivative
105                .index_axis_mut(Axis(0), i)
106                .assign(&phi.mapv(|p| p.eps1));
107            wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps1 = 0.0);
108        }
109        helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
110        Ok(())
111    }
112}
feos_dft/functional_contribution.rs

feos_dft/
functional_contribution.rs
