feos_dft/
functional.rs

1use crate::convolver::Convolver;
2use crate::functional_contribution::*;
3use crate::ideal_chain_contribution::IdealChainContribution;
4use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape};
5use feos_core::{EquationOfState, FeosResult, Residual, ResidualDyn, StateHD};
6use nalgebra::{DVector, dvector};
7use ndarray::*;
8use num_dual::*;
9use petgraph::Directed;
10use petgraph::graph::{Graph, UnGraph};
11use petgraph::visit::EdgeRef;
12use std::borrow::Cow;
13use std::ops::{Deref, MulAssign};
14
15impl<I: Clone, F: HelmholtzEnergyFunctionalDyn> HelmholtzEnergyFunctionalDyn
16    for EquationOfState<Vec<I>, F>
17{
18    type Contribution<'a>
19        = F::Contribution<'a>
20    where
21        Self: 'a;
22
23    fn contributions<'a>(&'a self) -> impl Iterator<Item = Self::Contribution<'a>> {
24        self.residual.contributions()
25    }
26
27    fn molecule_shape(&self) -> MoleculeShape<'_> {
28        self.residual.molecule_shape()
29    }
30
31    fn bond_lengths<N: DualNum<f64> + Copy>(&self, temperature: N) -> UnGraph<(), N> {
32        self.residual.bond_lengths(temperature)
33    }
34}
35
36/// Different representations for molecules within DFT.
37pub enum MoleculeShape<'a> {
38    /// For spherical molecules, the number of components.
39    Spherical(usize),
40    /// For non-spherical molecules in a homosegmented approach, the
41    /// chain length parameter $m$.
42    NonSpherical(&'a DVector<f64>),
43    /// For non-spherical molecules in a heterosegmented approach,
44    /// the component index for every segment.
45    Heterosegmented(&'a DVector<usize>),
46}
47
48/// A general Helmholtz energy functional.
49pub trait HelmholtzEnergyFunctionalDyn: ResidualDyn {
50    type Contribution<'a>: FunctionalContribution
51    where
52        Self: 'a;
53
54    /// Return a slice of [FunctionalContribution]s.
55    fn contributions<'a>(&'a self) -> impl Iterator<Item = Self::Contribution<'a>>;
56
57    /// Return the shape of the molecules and the necessary specifications.
58    fn molecule_shape(&self) -> MoleculeShape<'_>;
59
60    /// Overwrite this, if the functional consists of heterosegmented chains.
61    fn bond_lengths<N: DualNum<f64> + Copy>(&self, _temperature: N) -> UnGraph<(), N> {
62        Graph::with_capacity(0, 0)
63    }
64}
65
66/// A general Helmholtz energy functional.
67pub trait HelmholtzEnergyFunctional: Residual {
68    type Contribution<'a>: FunctionalContribution
69    where
70        Self: 'a;
71
72    /// Return a slice of [FunctionalContribution]s.
73    fn contributions<'a>(&'a self) -> impl Iterator<Item = Self::Contribution<'a>>;
74
75    /// Return the shape of the molecules and the necessary specifications.
76    fn molecule_shape(&self) -> MoleculeShape<'_>;
77
78    /// Overwrite this, if the functional consists of heterosegmented chains.
79    fn bond_lengths<N: DualNum<f64> + Copy>(&self, _temperature: N) -> UnGraph<(), N> {
80        Graph::with_capacity(0, 0)
81    }
82
83    fn weight_functions(&self, temperature: f64) -> Vec<WeightFunctionInfo<f64>> {
84        self.contributions()
85            .map(|c| c.weight_functions(temperature))
86            .collect()
87    }
88
89    fn m(&self) -> Cow<'_, [f64]> {
90        match self.molecule_shape() {
91            MoleculeShape::Spherical(n) => Cow::Owned(vec![1.0; n]),
92            MoleculeShape::NonSpherical(m) => Cow::Borrowed(m.as_slice()),
93            MoleculeShape::Heterosegmented(component_index) => {
94                Cow::Owned(vec![1.0; component_index.len()])
95            }
96        }
97    }
98
99    fn component_index(&self) -> Cow<'_, [usize]> {
100        match self.molecule_shape() {
101            MoleculeShape::Spherical(n) => Cow::Owned((0..n).collect()),
102            MoleculeShape::NonSpherical(m) => Cow::Owned((0..m.len()).collect()),
103            MoleculeShape::Heterosegmented(component_index) => {
104                Cow::Borrowed(component_index.as_slice())
105            }
106        }
107    }
108
109    fn ideal_chain_contribution(&self) -> IdealChainContribution {
110        IdealChainContribution::new(&self.component_index(), &self.m())
111    }
112
113    /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{\beta F}}{\delta\rho_i(\mathbf{r})}$.
114    #[expect(clippy::type_complexity)]
115    fn functional_derivative<D, N: DualNum<f64> + Copy>(
116        &self,
117        temperature: N,
118        density: &Array<N, D::Larger>,
119        convolver: &dyn Convolver<N, D>,
120    ) -> FeosResult<(Array<N, D>, Array<N, D::Larger>)>
121    where
122        D: Dimension,
123        D::Larger: Dimension<Smaller = D>,
124    {
125        let weighted_densities = convolver.weighted_densities(density);
126        let contributions = self.contributions();
127        let mut partial_derivatives = Vec::new();
128        let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
129        for (c, wd) in contributions.into_iter().zip(weighted_densities) {
130            let nwd = wd.shape()[0];
131            let ngrid = wd.len() / nwd;
132            let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
133            let mut pd = Array::zeros(wd.raw_dim());
134            c.first_partial_derivatives(
135                temperature,
136                wd.into_shape_with_order((nwd, ngrid)).unwrap(),
137                phi.view_mut().into_shape_with_order(ngrid).unwrap(),
138                pd.view_mut().into_shape_with_order((nwd, ngrid)).unwrap(),
139            )?;
140            partial_derivatives.push(pd);
141            helmholtz_energy_density += &phi;
142        }
143        Ok((
144            helmholtz_energy_density,
145            convolver.functional_derivative(&partial_derivatives),
146        ))
147    }
148
149    /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
150    fn bond_integrals<D, N: DualNum<f64> + Copy>(
151        &self,
152        temperature: N,
153        exponential: &Array<N, D::Larger>,
154        convolver: &dyn Convolver<N, D>,
155    ) -> Array<N, D::Larger>
156    where
157        D: Dimension,
158        D::Larger: Dimension<Smaller = D>,
159    {
160        // calculate weight functions
161        let bond_lengths = self.bond_lengths(temperature).into_edge_type();
162        let mut bond_weight_functions = bond_lengths.map(
163            |_, _| (),
164            |_, &l| WeightFunction::new_scaled(dvector![l], WeightFunctionShape::Delta),
165        );
166        for n in bond_lengths.node_indices() {
167            for e in bond_lengths.edges(n) {
168                bond_weight_functions.add_edge(
169                    e.target(),
170                    e.source(),
171                    WeightFunction::new_scaled(dvector![*e.weight()], WeightFunctionShape::Delta),
172                );
173            }
174        }
175
176        let mut i_graph: Graph<_, Option<Array<N, D>>, Directed> =
177            bond_weight_functions.map(|_, _| (), |_, _| None);
178
179        let bonds = i_graph.edge_count();
180        let mut calc = 0;
181
182        // go through the whole graph until every bond has been calculated
183        while calc < bonds {
184            let mut edge_id = None;
185            let mut i1 = None;
186
187            // find the first bond that can be calculated
188            'nodes: for node in i_graph.node_indices() {
189                for edge in i_graph.edges(node) {
190                    // skip already calculated bonds
191                    if edge.weight().is_some() {
192                        continue;
193                    }
194
195                    // if all bonds from the neighboring segment are calculated calculate the bond
196                    let edges = i_graph
197                        .edges(edge.target())
198                        .filter(|e| e.target().index() != node.index());
199                    if edges.clone().all(|e| e.weight().is_some()) {
200                        edge_id = Some(edge.id());
201                        let i0 = edges.fold(
202                            exponential
203                                .index_axis(Axis(0), edge.target().index())
204                                .to_owned(),
205                            |acc: Array<N, D>, e| acc * e.weight().as_ref().unwrap(),
206                        );
207                        i1 = Some(convolver.convolve(i0, &bond_weight_functions[edge.id()]));
208                        break 'nodes;
209                    }
210                }
211            }
212            if let Some(edge_id) = edge_id {
213                i_graph[edge_id] = i1;
214                calc += 1;
215            } else {
216                panic!("Cycle in molecular structure detected!")
217            }
218        }
219
220        let mut i = Array::ones(exponential.raw_dim());
221        for node in i_graph.node_indices() {
222            for edge in i_graph.edges(node) {
223                i.index_axis_mut(Axis(0), node.index())
224                    .mul_assign(edge.weight().as_ref().unwrap());
225            }
226        }
227
228        i
229    }
230
231    fn evaluate_bulk<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> Vec<(&'static str, D)> {
232        let mut res: Vec<_> = self
233            .contributions()
234            .map(|c| (c.name(), c.bulk_helmholtz_energy_density(state)))
235            .collect();
236        res.push((
237            "Ideal chain",
238            self.ideal_chain_contribution()
239                .bulk_helmholtz_energy_density(&state.partial_density),
240        ));
241        res
242    }
243}
244
245impl<C: Deref<Target = F> + Clone, F: HelmholtzEnergyFunctionalDyn + ResidualDyn + 'static>
246    HelmholtzEnergyFunctional for C
247{
248    type Contribution<'a>
249        = F::Contribution<'a>
250    where
251        Self: 'a;
252
253    fn contributions<'a>(&'a self) -> impl Iterator<Item = Self::Contribution<'a>> {
254        F::contributions(self.deref())
255    }
256
257    fn molecule_shape(&self) -> MoleculeShape<'_> {
258        F::molecule_shape(self.deref())
259    }
260
261    fn bond_lengths<N: DualNum<f64> + Copy>(&self, temperature: N) -> UnGraph<(), N> {
262        F::bond_lengths(self.deref(), temperature)
263    }
264}
feos_dft/functional.rs

feos_dft/
functional.rs
