feos_dft/
functional_contribution.rs

1use crate::weight_functions::WeightFunctionInfo;
2use feos_core::{FeosResult, StateHD};
3use ndarray::RemoveAxis;
4use ndarray::prelude::*;
5use num_dual::*;
6use num_traits::Zero;
7
8/// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers.
9pub trait FunctionalContribution: Sync + Send {
10    /// Return the name of the contribution.
11    fn name(&self) -> &'static str;
12
13    /// Return the weight functions required in this contribution.
14    fn weight_functions<N: DualNum<f64> + Copy>(&self, temperature: N) -> WeightFunctionInfo<N>;
15
16    /// Overwrite this if the weight functions in pDGT are different than for DFT.
17    fn weight_functions_pdgt<N: DualNum<f64> + Copy>(
18        &self,
19        temperature: N,
20    ) -> WeightFunctionInfo<N> {
21        self.weight_functions(temperature)
22    }
23
24    /// Return the Helmholtz energy density for the given temperature and weighted densities.
25    fn helmholtz_energy_density<N: DualNum<f64> + Copy>(
26        &self,
27        temperature: N,
28        weighted_densities: ArrayView2<N>,
29    ) -> FeosResult<Array1<N>>;
30
31    fn bulk_helmholtz_energy_density<N: DualNum<f64> + Copy>(&self, state: &StateHD<N>) -> N {
32        // calculate weight functions
33        let weight_functions = self.weight_functions(state.temperature);
34
35        // calculate segment density
36        let density: Array1<_> = weight_functions
37            .component_index
38            .iter()
39            .map(|&c| state.partial_density[c])
40            .collect();
41
42        // calculate weighted density and Helmholtz energy
43        let weight_constants = weight_functions.weight_constants(Zero::zero(), 0);
44        let weighted_densities = weight_constants.dot(&density).insert_axis(Axis(1));
45        self.helmholtz_energy_density(state.temperature, weighted_densities.view())
46            .unwrap()[0]
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    ) -> FeosResult<()> {
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    ) -> FeosResult<()> {
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
