feos-org · prehner · Sep 2, 2026 · Jul 26, 2026
diff --git a/crates/feos-dft/Cargo.toml b/crates/feos-dft/Cargo.toml
@@ -24,6 +24,7 @@ rustfft = { workspace = true }
 num-traits = { workspace = true }
 libm = { workspace = true }
 petgraph = { workspace = true }
+indexmap = { workspace = true }
 
 feos-core = { workspace = true }
 

diff --git a/crates/feos-dft/src/adsorption/mod.rs b/crates/feos-dft/src/adsorption/mod.rs
@@ -394,13 +394,13 @@ where
         })
     }
 
-    pub fn partial_molar_enthalpy_of_adsorption(&self) -> MolarEnergy<DMatrix<f64>> {
+    pub fn partial_molar_enthalpy_of_adsorption(&mut self) -> MolarEnergy<DMatrix<f64>> {
         let h_ads: Vec<_> = self
             .profiles
-            .iter()
+            .iter_mut()
             .map(|p| {
                 match p
-                    .as_ref()
+                    .as_mut()
                     .ok()
                     .and_then(|p| p.partial_molar_enthalpy_of_adsorption().ok())
                 {
@@ -414,10 +414,10 @@ where
         MolarEnergy::from_fn(self.components, self.profiles.len(), |j, i| h_ads[i].get(j))
     }
 
-    pub fn enthalpy_of_adsorption(&self) -> MolarEnergy<Array1<f64>> {
+    pub fn enthalpy_of_adsorption(&mut self) -> MolarEnergy<Array1<f64>> {
         MolarEnergy::from_shape_fn(self.profiles.len(), |i| {
             match self.profiles[i]
-                .as_ref()
+                .as_mut()
                 .ok()
                 .and_then(|p| p.enthalpy_of_adsorption().ok())
             {

diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs
@@ -123,7 +123,7 @@ where
         density: Option<&Density<Array<f64, D::Larger>>>,
         specification: PoreSpecification,
     ) -> Self {
-        let mut profile = DFTProfile::new(grid, bulk, Some(external_potential), density, Some(1));
+        let mut profile = DFTProfile::new(grid, bulk, Some(external_potential), density);
 
         // fix the number of particles
         match specification {
@@ -169,17 +169,19 @@ where
         self
     }
 
-    pub fn partial_molar_enthalpy_of_adsorption(&self) -> FeosResult<MolarEnergy<DVector<f64>>> {
+    pub fn partial_molar_enthalpy_of_adsorption(
+        &mut self,
+    ) -> FeosResult<MolarEnergy<DVector<f64>>> {
         let a = self.profile.dn_dmu()?;
         let a_unit = a.get2(0, 0);
-        let b = -self.profile.temperature * self.profile.dn_dt()?;
+        let b = -self.profile.bulk.temperature * self.profile.dn_dt()?;
         let b_unit = b.get(0);
 
         let h_ads = LU::new((a / a_unit).into_value())?.solve(&(b / b_unit).into_value());
         Ok(&h_ads * b_unit / a_unit)
     }
 
-    pub fn enthalpy_of_adsorption(&self) -> FeosResult<MolarEnergy> {
+    pub fn enthalpy_of_adsorption(&mut self) -> FeosResult<MolarEnergy> {
         Ok(self
             .partial_molar_enthalpy_of_adsorption()?
             .dot(&Dimensionless::new(self.profile.bulk.molefracs.clone())))
@@ -195,16 +197,15 @@ where
             )
         };
         let pot = (self.profile.external_potential.mapv(N::from)
-            * self.profile.temperature.to_reduced())
+            * self.profile.bulk.temperature.to_reduced())
         .mapv(|v| v / temperature);
         let exp_pot = pot.mapv(|v| (-v).exp());
         let functional_contributions = self.profile.bulk.eos.contributions();
         let weight_functions: Vec<WeightFunctionInfo<N>> = functional_contributions
             .into_iter()
             .map(|c| c.weight_functions(temperature))
             .collect();
-        let convolver =
-            ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions, self.profile.lanczos);
+        let convolver = ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions);
         let bonds = self
             .profile
             .bulk
@@ -214,38 +215,35 @@ where
     }
 
     pub fn henry_coefficients(&self) -> HenryCoefficient<DVector<f64>> {
-        let t = self.profile.temperature.to_reduced();
-        Volume::from_reduced(self._henry_coefficients(t)) / (RGAS * self.profile.temperature)
+        let t = self.profile.bulk.temperature.to_reduced();
+        Volume::from_reduced(self._henry_coefficients(t)) / (RGAS * self.profile.bulk.temperature)
     }
 
     pub fn ideal_gas_enthalpy_of_adsorption(&self) -> MolarEnergy<DVector<f64>> {
-        let t = Dual64::from(self.profile.temperature.to_reduced()).derivative();
+        let t = Dual64::from(self.profile.bulk.temperature.to_reduced()).derivative();
         let h_dual = self._henry_coefficients(t);
         let h = h_dual.map(|h| h.re);
         let dh = h_dual.map(|h| h.eps);
-        let t = self.profile.temperature.to_reduced();
-        RGAS * self.profile.temperature
+        let t = self.profile.bulk.temperature.to_reduced();
+        RGAS * self.profile.bulk.temperature
             * Dimensionless::from_reduced((&h - t * dh).component_div(&h))
     }
 
     pub fn into_dyn(self) -> PoreProfile<IxDyn, F> {
         // initialize convolver
         let t = self.profile.bulk.temperature.to_reduced();
         let weight_functions = self.profile.bulk.eos.weight_functions(t);
-        let convolver =
-            ConvolverFFT::plan(&self.profile.grid, &weight_functions, self.profile.lanczos);
+        let convolver = ConvolverFFT::plan(&self.profile.grid, &weight_functions);
 
         PoreProfile {
             profile: DFTProfile {
                 grid: self.profile.grid,
                 convolver,
-                temperature: self.profile.temperature,
                 density: self.profile.density.into_dyn(),
                 specification: self.profile.specification,
                 external_potential: self.profile.external_potential.into_dyn(),
                 bulk: self.profile.bulk,
                 solver_log: self.profile.solver_log,
-                lanczos: self.profile.lanczos,
             },
             grand_potential: self.grand_potential,
             interfacial_tension: self.interfacial_tension,

diff --git a/crates/feos-dft/src/convolver/mod.rs b/crates/feos-dft/src/convolver/mod.rs
@@ -93,8 +93,9 @@ where
     pub fn plan(
         grid: &Grid,
         weight_functions: &[WeightFunctionInfo<T>],
-        lanczos: Option<i32>,
     ) -> Arc<dyn Convolver<T, D>> {
+        // For consistency, we always use an exponent of 1.
+        let lanczos = Some(1);
         match grid {
             Grid::Bulk => PeriodicConvolver::new_0d(weight_functions),
             Grid::Polar(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos),

diff --git a/crates/feos-dft/src/functional.rs b/crates/feos-dft/src/functional.rs
@@ -1,16 +1,16 @@
 use crate::convolver::Convolver;
-use crate::functional_contribution::*;
 use crate::ideal_chain_contribution::IdealChainContribution;
 use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape};
-use feos_core::{EquationOfState, FeosResult, Residual, ResidualDyn, StateHD};
+use crate::{DFTSolverLog, functional_contribution::*};
+use feos_core::{EquationOfState, FeosError, FeosResult, Residual, ResidualDyn, StateHD};
 use nalgebra::{DVector, dvector};
 use ndarray::*;
 use num_dual::*;
 use petgraph::Directed;
 use petgraph::graph::{Graph, UnGraph};
 use petgraph::visit::EdgeRef;
 use std::borrow::Cow;
-use std::ops::{Deref, MulAssign};
+use std::ops::{AddAssign, Deref, MulAssign};
 
 impl<I: Clone, F: HelmholtzEnergyFunctionalDyn> HelmholtzEnergyFunctionalDyn
     for EquationOfState<Vec<I>, F>
@@ -112,38 +112,87 @@ pub trait HelmholtzEnergyFunctional: Residual {
 
     /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{\beta F}}{\delta\rho_i(\mathbf{r})}$.
     #[expect(clippy::type_complexity)]
-    fn functional_derivative<D, N: DualNum<Primitive = f64> + Copy>(
+    fn functional_derivative<D, N: DualNumCopy<Primitive = f64>>(
         &self,
         temperature: N,
         density: &Array<N, D::Larger>,
         convolver: &dyn Convolver<N, D>,
+        solver_log: &mut DFTSolverLog,
     ) -> FeosResult<(Array<N, D>, Array<N, D::Larger>)>
     where
         D: Dimension,
         D::Larger: Dimension<Smaller = D>,
     {
-        let weighted_densities = convolver.weighted_densities(density);
+        // calculate weighted densities
+        let weighted_densities = solver_log.time_function("weighted densities", || {
+            convolver.weighted_densities(density)
+        });
+
+        // calculate partial derivatives
         let contributions = self.contributions();
         let mut partial_derivatives = Vec::new();
         let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
-        for (c, wd) in contributions.into_iter().zip(weighted_densities) {
+        solver_log.time_function("partial derivatives", || {
+            for (c, wd) in contributions.into_iter().zip(weighted_densities) {
+                let nwd = wd.shape()[0];
+                let ngrid = wd.len() / nwd;
+                let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
+                let mut pd = Array::zeros(wd.raw_dim());
+                c.first_partial_derivatives(
+                    temperature,
+                    wd.into_shape_with_order((nwd, ngrid)).unwrap(),
+                    phi.view_mut().into_shape_with_order(ngrid).unwrap(),
+                    pd.view_mut().into_shape_with_order((nwd, ngrid)).unwrap(),
+                )?;
+                partial_derivatives.push(pd);
+                helmholtz_energy_density += &phi;
+            }
+            Ok::<_, FeosError>(())
+        })?;
+
+        // calculate functional derivative
+        let functional_derivative = solver_log.time_function("functional derivative", || {
+            convolver.functional_derivative(&partial_derivatives)
+        });
+
+        Ok((helmholtz_energy_density, functional_derivative))
+    }
+
+    fn second_partial_derivatives<D>(
+        &self,
+        temperature: f64,
+        density: &Array<f64, D::Larger>,
+        convolver: &dyn Convolver<f64, D>,
+    ) -> FeosResult<Vec<Array<f64, <D::Larger as Dimension>::Larger>>>
+    where
+        D: Dimension,
+        D::Larger: Dimension<Smaller = D>,
+    {
+        let contributions = self.contributions();
+
+        let weighted_densities = convolver.weighted_densities(density);
+
+        let mut second_partial_derivatives = Vec::new();
+        for (c, wd) in contributions.into_iter().zip(&weighted_densities) {
             let nwd = wd.shape()[0];
             let ngrid = wd.len() / nwd;
             let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
             let mut pd = Array::zeros(wd.raw_dim());
-            c.first_partial_derivatives(
+            let dim = wd.shape();
+            let dim: Vec<_> = std::iter::once(&nwd).chain(dim).cloned().collect();
+            let mut pd2 = Array::zeros(dim).into_dimensionality().unwrap();
+            c.second_partial_derivatives(
                 temperature,
-                wd.into_shape_with_order((nwd, ngrid)).unwrap(),
+                wd.view().into_shape_with_order((nwd, ngrid)).unwrap(),
                 phi.view_mut().into_shape_with_order(ngrid).unwrap(),
                 pd.view_mut().into_shape_with_order((nwd, ngrid)).unwrap(),
+                pd2.view_mut()
+                    .into_shape_with_order((nwd, nwd, ngrid))
+                    .unwrap(),
             )?;
-            partial_derivatives.push(pd);
-            helmholtz_energy_density += &phi;
+            second_partial_derivatives.push(pd2);
         }
-        Ok((
-            helmholtz_energy_density,
-            convolver.functional_derivative(&partial_derivatives),
-        ))
+        Ok(second_partial_derivatives)
     }
 
     /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
@@ -228,6 +277,106 @@ pub trait HelmholtzEnergyFunctional: Residual {
         i
     }
 
+    fn delta_bond_integrals<D>(
+        &self,
+        temperature: f64,
+        exponential: &Array<f64, D::Larger>,
+        delta_functional_derivative: &Array<f64, D::Larger>,
+        convolver: &dyn Convolver<f64, D>,
+    ) -> Array<f64, D::Larger>
+    where
+        D: Dimension,
+        D::Larger: Dimension<Smaller = D>,
+    {
+        // calculate weight functions
+        let bond_lengths = self.bond_lengths(temperature).into_edge_type();
+        let mut bond_weight_functions = bond_lengths.map(
+            |_, _| (),
+            |_, &l| WeightFunction::new_scaled(dvector![l], WeightFunctionShape::Delta),
+        );
+        for n in bond_lengths.node_indices() {
+            for e in bond_lengths.edges(n) {
+                bond_weight_functions.add_edge(
+                    e.target(),
+                    e.source(),
+                    WeightFunction::new_scaled(dvector![*e.weight()], WeightFunctionShape::Delta),
+                );
+            }
+        }
+
+        let mut i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
+            bond_weight_functions.map(|_, _| (), |_, _| None);
+        let mut delta_i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
+            bond_weight_functions.map(|_, _| (), |_, _| None);
+
+        let bonds = i_graph.edge_count();
+        let mut calc = 0;
+
+        // go through the whole graph until every bond has been calculated
+        while calc < bonds {
+            let mut edge_id = None;
+            let mut i1 = None;
+            let mut delta_i1 = None;
+
+            // find the first bond that can be calculated
+            'nodes: for node in i_graph.node_indices() {
+                for edge in i_graph.edges(node) {
+                    // skip already calculated bonds
+                    if edge.weight().is_some() {
+                        continue;
+                    }
+
+                    // if all bonds from the neighboring segment are calculated calculate the bond
+                    let edges = i_graph
+                        .edges(edge.target())
+                        .filter(|e| e.target().index() != node.index());
+                    let delta_edges = delta_i_graph
+                        .edges(edge.target())
+                        .filter(|e| e.target().index() != node.index());
+                    if edges.clone().all(|e| e.weight().is_some()) {
+                        edge_id = Some(edge.id());
+                        let i0 = edges.fold(
+                            exponential
+                                .index_axis(Axis(0), edge.target().index())
+                                .to_owned(),
+                            |acc: Array<f64, _>, e| acc * e.weight().as_ref().unwrap(),
+                        );
+                        let delta_i0 = delta_edges.fold(
+                            -&delta_functional_derivative
+                                .index_axis(Axis(0), edge.target().index()),
+                            |acc: Array<f64, _>, delta_e| acc + delta_e.weight().as_ref().unwrap(),
+                        ) * &i0;
+                        i1 = Some(convolver.convolve(i0, &bond_weight_functions[edge.id()]));
+                        delta_i1 = Some(
+                            (convolver.convolve(delta_i0, &bond_weight_functions[edge.id()])
+                                / i1.as_ref().unwrap())
+                            .mapv(|x| if x.is_finite() { x } else { 0.0 }),
+                        );
+                        break 'nodes;
+                    }
+                }
+            }
+            if let Some(edge_id) = edge_id {
+                i_graph[edge_id] = i1;
+                delta_i_graph[edge_id] = delta_i1;
+                calc += 1;
+            } else {
+                panic!("Cycle in molecular structure detected!")
+            }
+        }
+
+        let mut delta_i = Array::zeros(exponential.raw_dim());
+        for node in delta_i_graph.node_indices() {
+            for edge in delta_i_graph.edges(node) {
+                delta_i
+                    .index_axis_mut(Axis(0), node.index())
+                    .add_assign(edge.weight().as_ref().unwrap());
+            }
+        }
+
+        delta_i
+    }
+
     fn evaluate_bulk<D: DualNum<Primitive = f64> + Copy>(
         &self,
         state: &StateHD<D>,