feos_dft/
functional.rs

1use crate::adsorption::FluidParameters;
2use crate::convolver::Convolver;
3use crate::functional_contribution::*;
4use crate::ideal_chain_contribution::IdealChainContribution;
5use crate::solvation::PairPotential;
6use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape};
7use feos_core::{Components, EosResult, EquationOfState, IdealGas, Molarweight, Residual, StateHD};
8use ndarray::*;
9use num_dual::*;
10use petgraph::graph::{Graph, UnGraph};
11use petgraph::visit::EdgeRef;
12use petgraph::Directed;
13use quantity::MolarWeight;
14use std::borrow::Cow;
15use std::ops::{Deref, MulAssign};
16use std::sync::Arc;
17
18impl<I: Components + Send + Sync, F: HelmholtzEnergyFunctional> HelmholtzEnergyFunctional
19    for EquationOfState<I, F>
20{
21    type Contribution = F::Contribution;
22
23    fn contributions(&self) -> Box<dyn Iterator<Item = Self::Contribution>> {
24        self.residual.contributions()
25    }
26
27    fn molecule_shape(&self) -> MoleculeShape {
28        self.residual.molecule_shape()
29    }
30
31    fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
32        self.residual.compute_max_density(moles)
33    }
34
35    fn bond_lengths<N: DualNum<f64> + Copy>(&self, temperature: N) -> UnGraph<(), N> {
36        self.residual.bond_lengths(temperature)
37    }
38}
39
40impl<I, F: PairPotential> PairPotential for EquationOfState<I, F> {
41    fn pair_potential(&self, i: usize, r: &Array1<f64>, temperature: f64) -> Array2<f64> {
42        self.residual.pair_potential(i, r, temperature)
43    }
44}
45
46impl<I: Components + Send + Sync, F: FluidParameters> FluidParameters for EquationOfState<I, F> {
47    fn epsilon_k_ff(&self) -> Array1<f64> {
48        self.residual.epsilon_k_ff()
49    }
50
51    fn sigma_ff(&self) -> &Array1<f64> {
52        self.residual.sigma_ff()
53    }
54}
55
56/// Wrapper struct for the [HelmholtzEnergyFunctional] trait.
57///
58/// Needed (for now) to generically implement the `Residual`
59/// trait for Helmholtz energy functionals.
60#[derive(Clone)]
61pub struct DFT<F>(pub F);
62
63impl<F> DFT<F> {
64    pub fn into<F2: From<F>>(self) -> DFT<F2> {
65        DFT(self.0.into())
66    }
67}
68
69impl<F> Deref for DFT<F> {
70    type Target = F;
71    fn deref(&self) -> &F {
72        &self.0
73    }
74}
75
76impl<F> DFT<F> {
77    pub fn ideal_gas<I>(self, ideal_gas: I) -> DFT<EquationOfState<I, F>> {
78        DFT(EquationOfState::new(Arc::new(ideal_gas), Arc::new(self.0)))
79    }
80}
81
82impl<F: HelmholtzEnergyFunctional> Components for DFT<F> {
83    fn components(&self) -> usize {
84        self.0.components()
85    }
86
87    fn subset(&self, component_list: &[usize]) -> Self {
88        Self(self.0.subset(component_list))
89    }
90}
91
92impl<F: HelmholtzEnergyFunctional> Residual for DFT<F> {
93    fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
94        self.0.compute_max_density(moles)
95    }
96
97    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
98        &self,
99        state: &StateHD<D>,
100    ) -> Vec<(String, D)> {
101        let mut res: Vec<(String, D)> = self
102            .0
103            .contributions()
104            .map(|c| (c.to_string(), c.helmholtz_energy(state)))
105            .collect();
106        res.push((
107            self.ideal_chain_contribution().to_string(),
108            self.ideal_chain_contribution().helmholtz_energy(state),
109        ));
110        res
111    }
112}
113
114impl<F: Molarweight> Molarweight for DFT<F> {
115    fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
116        self.0.molar_weight()
117    }
118}
119
120impl<F: HelmholtzEnergyFunctional + IdealGas> IdealGas for DFT<F> {
121    fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
122        self.0.ln_lambda3(temperature)
123    }
124
125    fn ideal_gas_model(&self) -> String {
126        self.0.ideal_gas_model()
127    }
128}
129
130/// Different representations for molecules within DFT.
131pub enum MoleculeShape<'a> {
132    /// For spherical molecules, the number of components.
133    Spherical(usize),
134    /// For non-spherical molecules in a homosegmented approach, the
135    /// chain length parameter $m$.
136    NonSpherical(&'a Array1<f64>),
137    /// For non-spherical molecules in a heterosegmented approach,
138    /// the component index for every segment.
139    Heterosegmented(&'a Array1<usize>),
140}
141
142/// A general Helmholtz energy functional.
143pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
144    type Contribution: FunctionalContribution;
145
146    /// Return a slice of [FunctionalContribution]s.
147    fn contributions(&self) -> Box<dyn Iterator<Item = Self::Contribution>>;
148
149    /// Return the shape of the molecules and the necessary specifications.
150    fn molecule_shape(&self) -> MoleculeShape;
151
152    /// Return the maximum density in Angstrom^-3.
153    ///
154    /// This value is used as an estimate for a liquid phase for phase
155    /// equilibria and other iterations. It is not explicitly meant to
156    /// be a mathematical limit for the density (if those exist in the
157    /// equation of state anyways).
158    fn compute_max_density(&self, moles: &Array1<f64>) -> f64;
159
160    /// Overwrite this, if the functional consists of heterosegmented chains.
161    fn bond_lengths<N: DualNum<f64> + Copy>(&self, _temperature: N) -> UnGraph<(), N> {
162        Graph::with_capacity(0, 0)
163    }
164
165    fn weight_functions(&self, temperature: f64) -> Vec<WeightFunctionInfo<f64>> {
166        self.contributions()
167            .map(|c| c.weight_functions(temperature))
168            .collect()
169    }
170
171    fn m(&self) -> Cow<Array1<f64>> {
172        match self.molecule_shape() {
173            MoleculeShape::Spherical(n) => Cow::Owned(Array1::ones(n)),
174            MoleculeShape::NonSpherical(m) => Cow::Borrowed(m),
175            MoleculeShape::Heterosegmented(component_index) => {
176                Cow::Owned(Array1::ones(component_index.len()))
177            }
178        }
179    }
180
181    fn component_index(&self) -> Cow<Array1<usize>> {
182        match self.molecule_shape() {
183            MoleculeShape::Spherical(n) => Cow::Owned(Array1::from_shape_fn(n, |i| i)),
184            MoleculeShape::NonSpherical(m) => Cow::Owned(Array1::from_shape_fn(m.len(), |i| i)),
185            MoleculeShape::Heterosegmented(component_index) => Cow::Borrowed(component_index),
186        }
187    }
188
189    fn ideal_chain_contribution(&self) -> IdealChainContribution {
190        IdealChainContribution::new(&self.component_index(), &self.m())
191    }
192
193    /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{\beta F}}{\delta\rho_i(\mathbf{r})}$.
194    #[expect(clippy::type_complexity)]
195    fn functional_derivative<D, N: DualNum<f64> + Copy + ScalarOperand>(
196        &self,
197        temperature: N,
198        density: &Array<N, D::Larger>,
199        convolver: &Arc<dyn Convolver<N, D>>,
200    ) -> EosResult<(Array<N, D>, Array<N, D::Larger>)>
201    where
202        D: Dimension,
203        D::Larger: Dimension<Smaller = D>,
204    {
205        let weighted_densities = convolver.weighted_densities(density);
206        let contributions = self.contributions();
207        let mut partial_derivatives = Vec::new();
208        let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
209        for (c, wd) in contributions.zip(weighted_densities) {
210            let nwd = wd.shape()[0];
211            let ngrid = wd.len() / nwd;
212            let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
213            let mut pd = Array::zeros(wd.raw_dim());
214            c.first_partial_derivatives(
215                temperature,
216                wd.into_shape_with_order((nwd, ngrid)).unwrap(),
217                phi.view_mut().into_shape_with_order(ngrid).unwrap(),
218                pd.view_mut().into_shape_with_order((nwd, ngrid)).unwrap(),
219            )?;
220            partial_derivatives.push(pd);
221            helmholtz_energy_density += &phi;
222        }
223        Ok((
224            helmholtz_energy_density,
225            convolver.functional_derivative(&partial_derivatives),
226        ))
227    }
228
229    /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
230    fn bond_integrals<D, N: DualNum<f64> + Copy>(
231        &self,
232        temperature: N,
233        exponential: &Array<N, D::Larger>,
234        convolver: &Arc<dyn Convolver<N, D>>,
235    ) -> Array<N, D::Larger>
236    where
237        D: Dimension,
238        D::Larger: Dimension<Smaller = D>,
239    {
240        // calculate weight functions
241        let bond_lengths = self.bond_lengths(temperature).into_edge_type();
242        let mut bond_weight_functions = bond_lengths.map(
243            |_, _| (),
244            |_, &l| WeightFunction::new_scaled(arr1(&[l]), WeightFunctionShape::Delta),
245        );
246        for n in bond_lengths.node_indices() {
247            for e in bond_lengths.edges(n) {
248                bond_weight_functions.add_edge(
249                    e.target(),
250                    e.source(),
251                    WeightFunction::new_scaled(arr1(&[*e.weight()]), WeightFunctionShape::Delta),
252                );
253            }
254        }
255
256        let mut i_graph: Graph<_, Option<Array<N, D>>, Directed> =
257            bond_weight_functions.map(|_, _| (), |_, _| None);
258
259        let bonds = i_graph.edge_count();
260        let mut calc = 0;
261
262        // go through the whole graph until every bond has been calculated
263        while calc < bonds {
264            let mut edge_id = None;
265            let mut i1 = None;
266
267            // find the first bond that can be calculated
268            'nodes: for node in i_graph.node_indices() {
269                for edge in i_graph.edges(node) {
270                    // skip already calculated bonds
271                    if edge.weight().is_some() {
272                        continue;
273                    }
274
275                    // if all bonds from the neighboring segment are calculated calculate the bond
276                    let edges = i_graph
277                        .edges(edge.target())
278                        .filter(|e| e.target().index() != node.index());
279                    if edges.clone().all(|e| e.weight().is_some()) {
280                        edge_id = Some(edge.id());
281                        let i0 = edges.fold(
282                            exponential
283                                .index_axis(Axis(0), edge.target().index())
284                                .to_owned(),
285                            |acc: Array<N, D>, e| acc * e.weight().as_ref().unwrap(),
286                        );
287                        i1 = Some(convolver.convolve(i0, &bond_weight_functions[edge.id()]));
288                        break 'nodes;
289                    }
290                }
291            }
292            if let Some(edge_id) = edge_id {
293                i_graph[edge_id] = i1;
294                calc += 1;
295            } else {
296                panic!("Cycle in molecular structure detected!")
297            }
298        }
299
300        let mut i = Array::ones(exponential.raw_dim());
301        for node in i_graph.node_indices() {
302            for edge in i_graph.edges(node) {
303                i.index_axis_mut(Axis(0), node.index())
304                    .mul_assign(edge.weight().as_ref().unwrap());
305            }
306        }
307
308        i
309    }
310}
feos_dft/functional.rs

feos_dft/
functional.rs
