Repository navigation
Expand file tree
/
Copy pathfunctional_contribution.rs
More file actions
112 lines (101 loc) · 4.38 KB
/
Copy pathfunctional_contribution.rs
File metadata and controls
112 lines (101 loc) · 4.38 KB
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
use crate::weight_functions::WeightFunctionInfo;
use feos_core::{FeosResult, StateHD};
use ndarray::RemoveAxis;
use ndarray::prelude::*;
use num_dual::*;
use num_traits::Zero;
/// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers.
pub trait FunctionalContribution: Sync + Send {
/// Return the name of the contribution.
fn name(&self) -> &'static str;
/// Return the weight functions required in this contribution.
fn weight_functions<N: DualNum<f64> + Copy>(&self, temperature: N) -> WeightFunctionInfo<N>;
/// Overwrite this if the weight functions in pDGT are different than for DFT.
fn weight_functions_pdgt<N: DualNum<f64> + Copy>(
&self,
temperature: N,
) -> WeightFunctionInfo<N> {
self.weight_functions(temperature)
}
/// Return the Helmholtz energy density for the given temperature and weighted densities.
fn helmholtz_energy_density<N: DualNum<f64> + Copy>(
&self,
temperature: N,
weighted_densities: ArrayView2<N>,
) -> FeosResult<Array1<N>>;
fn bulk_helmholtz_energy_density<N: DualNum<f64> + Copy>(&self, state: &StateHD<N>) -> N {
// calculate weight functions
let weight_functions = self.weight_functions(state.temperature);
// calculate segment density
let density: Array1<_> = weight_functions
.component_index
.iter()
.map(|&c| state.partial_density[c])
.collect();
// calculate weighted density and Helmholtz energy
let weight_constants = weight_functions.weight_constants(Zero::zero(), 0);
let weighted_densities = weight_constants.dot(&density).insert_axis(Axis(1));
self.helmholtz_energy_density(state.temperature, weighted_densities.view())
.unwrap()[0]
}
fn first_partial_derivatives<N: DualNum<f64> + Copy>(
&self,
temperature: N,
weighted_densities: Array2<N>,
mut helmholtz_energy_density: ArrayViewMut1<N>,
mut first_partial_derivative: ArrayViewMut2<N>,
) -> FeosResult<()> {
let mut wd = weighted_densities.mapv(Dual::from_re);
let t = Dual::from_re(temperature);
let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
for i in 0..wd.shape()[0] {
wd.index_axis_mut(Axis(0), i)
.map_inplace(|x| x.eps = N::one());
phi = self.helmholtz_energy_density(t, wd.view())?;
first_partial_derivative
.index_axis_mut(Axis(0), i)
.assign(&phi.mapv(|p| p.eps));
wd.index_axis_mut(Axis(0), i)
.map_inplace(|x| x.eps = N::zero());
}
helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
Ok(())
}
fn second_partial_derivatives(
&self,
temperature: f64,
weighted_densities: ArrayView2<f64>,
mut helmholtz_energy_density: ArrayViewMut1<f64>,
mut first_partial_derivative: ArrayViewMut2<f64>,
mut second_partial_derivative: ArrayViewMut3<f64>,
) -> FeosResult<()> {
let mut wd = weighted_densities.mapv(HyperDual64::from);
let t = HyperDual64::from(temperature);
let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
for i in 0..wd.shape()[0] {
wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps1 = 1.0);
for j in 0..=i {
wd.index_axis_mut(Axis(0), j).map_inplace(|x| x.eps2 = 1.0);
phi = self.helmholtz_energy_density(t, wd.view())?;
let p = phi.mapv(|p| p.eps1eps2);
second_partial_derivative
.index_axis_mut(Axis(0), i)
.index_axis_mut(Axis(0), j)
.assign(&p);
if i != j {
second_partial_derivative
.index_axis_mut(Axis(0), j)
.index_axis_mut(Axis(0), i)
.assign(&p);
}
wd.index_axis_mut(Axis(0), j).map_inplace(|x| x.eps2 = 0.0);
}
first_partial_derivative
.index_axis_mut(Axis(0), i)
.assign(&phi.mapv(|p| p.eps1));
wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps1 = 0.0);
}
helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
Ok(())
}
}