feos-org · prehner · Aug 20, 2026 · Aug 14, 2026 · Jul 26, 2026
diff --git a/crates/feos-dft/Cargo.toml b/crates/feos-dft/Cargo.toml
@@ -16,7 +16,7 @@ features = ["rayon"]
 
 [dependencies]
 quantity = { workspace = true, features = ["ndarray"] }
-num-dual = { workspace = true }
+num-dual = { workspace = true, features = ["ndarray"] }
 ndarray = { workspace = true }
 nalgebra = { workspace = true }
 rustdct = { workspace = true }

diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs
@@ -162,6 +162,7 @@ where
     }
 
     pub fn update_bulk(mut self, bulk: &State<F>) -> Self {
+        self.profile.specification = DFTSpecification::from_state(bulk);
         self.profile.bulk = bulk.clone();
         self.grand_potential = None;
         self.interfacial_tension = None;

diff --git a/crates/feos-dft/src/convolver/mod.rs b/crates/feos-dft/src/convolver/mod.rs
@@ -1,10 +1,8 @@
 use crate::geometry::{Axis, Geometry, Grid};
 use crate::weight_functions::*;
-use ndarray::linalg::Dot;
 use ndarray::prelude::*;
 use ndarray::{Axis as Axis_nd, RemoveAxis, Slice};
 use num_dual::*;
-use num_traits::Zero;
 use rustdct::DctNum;
 use std::ops::{AddAssign, MulAssign, SubAssign};
 use std::sync::Arc;
@@ -36,46 +34,6 @@ pub trait Convolver<T, D: Dimension>: Send + Sync {
     ) -> Array<T, D::Larger>;
 }
 
-pub(crate) struct BulkConvolver<T> {
-    weight_constants: Vec<Array2<T>>,
-}
-
-impl<T: DualNum<Primitive = f64> + Copy + Send + Sync> BulkConvolver<T> {
-    #[expect(clippy::new_ret_no_self)]
-    pub(crate) fn new(weight_functions: Vec<WeightFunctionInfo<T>>) -> Arc<dyn Convolver<T, Ix0>> {
-        let weight_constants = weight_functions
-            .into_iter()
-            .map(|w| w.weight_constants(Zero::zero(), 0))
-            .collect();
-        Arc::new(Self { weight_constants })
-    }
-}
-
-impl<T: DualNum<Primitive = f64> + Copy + Send + Sync> Convolver<T, Ix0> for BulkConvolver<T>
-where
-    Array2<T>: Dot<Array1<T>, Output = Array1<T>>,
-{
-    fn convolve(&self, _: Array0<T>, _: &WeightFunction<T>) -> Array0<T> {
-        unreachable!()
-    }
-
-    fn weighted_densities(&self, density: &Array1<T>) -> Vec<Array1<T>> {
-        self.weight_constants
-            .iter()
-            .map(|w| w.dot(density))
-            .collect()
-    }
-
-    fn functional_derivative(&self, partial_derivatives: &[Array1<T>]) -> Array1<T> {
-        self.weight_constants
-            .iter()
-            .zip(partial_derivatives.iter())
-            .map(|(w, pd)| pd.dot(w))
-            .reduce(|a, b| a + b)
-            .unwrap()
-    }
-}
-
 /// Base structure to hold either information about the weight function through
 /// `WeightFunctionInfo` or the weight functions themselves via
 /// `FFTWeightFunctions`.
@@ -125,11 +83,10 @@ pub struct ConvolverFFT<T, D: Dimension> {
     cartesian_transforms: Vec<CartesianTransform<T>>,
 }
 
-impl<T, D: Dimension + RemoveAxis + 'static> ConvolverFFT<T, D>
+impl<T, D: Dimension + 'static> ConvolverFFT<T, D>
 where
     T: DctNum + DualNum<Primitive = f64>,
     D::Larger: Dimension<Smaller = D>,
-    D::Smaller: Dimension<Larger = D>,
     <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
     /// Create the appropriate FFT convolver for the given grid.
@@ -139,9 +96,11 @@ where
         lanczos: Option<i32>,
     ) -> Arc<dyn Convolver<T, D>> {
         match grid {
+            Grid::Bulk => PeriodicConvolver::new_0d(weight_functions),
             Grid::Polar(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos),
             Grid::Spherical(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos),
             Grid::Cartesian1(z) => Self::new(Some(z), &[], weight_functions, lanczos),
+            Grid::Periodical1(z) => PeriodicConvolver::new_1d(z, weight_functions, lanczos),
             Grid::Cylindrical { r, z } => {
                 CurvilinearConvolver::new(r, &[z], weight_functions, lanczos)
             }
@@ -553,11 +512,10 @@ struct CurvilinearConvolver<T, D> {
     convolver_boundary: Arc<dyn Convolver<T, D>>,
 }
 
-impl<T, D: Dimension + RemoveAxis + 'static> CurvilinearConvolver<T, D>
+impl<T, D: Dimension + 'static> CurvilinearConvolver<T, D>
 where
     T: DctNum + DualNum<Primitive = f64>,
     D::Larger: Dimension<Smaller = D>,
-    D::Smaller: Dimension<Larger = D>,
     <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
     #[expect(clippy::new_ret_no_self)]
@@ -574,10 +532,9 @@ where
     }
 }
 
-impl<T, D: Dimension + RemoveAxis> Convolver<T, D> for CurvilinearConvolver<T, D>
+impl<T, D: Dimension> Convolver<T, D> for CurvilinearConvolver<T, D>
 where
     T: DctNum + DualNum<Primitive = f64>,
-    D::Smaller: Dimension<Larger = D>,
     D::Larger: Dimension<Smaller = D>,
 {
     fn convolve(
@@ -587,24 +544,22 @@ where
     ) -> Array<T, D> {
         // subtract boundary profile from full profile
         let profile_boundary = profile
-            .index_axis(Axis(0), profile.shape()[0] - 1)
+            .slice_axis(Axis(0), Slice::from(profile.shape()[0] - 1..))
             .into_owned();
-        for mut lane in profile.outer_iter_mut() {
+        for mut lane in profile.axis_chunks_iter_mut(Axis(0), 1) {
             lane.sub_assign(&profile_boundary);
         }
 
         // convolve full profile
         let mut result = self.convolver.convolve(profile, weight_function);
 
         // convolve boundary profile
-        let profile_boundary = profile_boundary.insert_axis(Axis(0));
         let result_boundary = self
             .convolver_boundary
             .convolve(profile_boundary, weight_function);
 
         // Add boundary result back to full result
-        let result_boundary = result_boundary.index_axis(Axis(0), 0);
-        for mut lane in result.outer_iter_mut() {
+        for mut lane in result.axis_chunks_iter_mut(Axis(0), 1) {
             lane.add_assign(&result_boundary);
         }
         result

diff --git a/crates/feos-dft/src/convolver/periodic_convolver.rs b/crates/feos-dft/src/convolver/periodic_convolver.rs
@@ -30,6 +30,18 @@ where
     D::Larger: Dimension<Smaller = D>,
     <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
+    pub fn new_0d(weight_functions: &[WeightFunctionInfo<T>]) -> Arc<dyn Convolver<T, D>> {
+        Self::new(&[], |_| {}, weight_functions, None)
+    }
+
+    pub fn new_1d(
+        axis: &Axis,
+        weight_functions: &[WeightFunctionInfo<T>],
+        lanczos: Option<i32>,
+    ) -> Arc<dyn Convolver<T, D>> {
+        Self::new(&[axis], |_| {}, weight_functions, lanczos)
+    }
+
     pub fn new_2d(
         axes: &[&Axis],
         angle: Angle,
@@ -200,7 +212,7 @@ impl<T: FftNum, D: Dimension> PeriodicConvolver<T, D> {
     }
 
     fn forward_transform<D2: Dimension>(&self, f: ArrayView<T, D2>) -> Array<Complex<T>, D2> {
-        let offset = D2::NDIM.unwrap() - D::NDIM.unwrap();
+        let offset = f.ndim() - self.k_abs.ndim();
         let mut result = f.mapv(Complex::from);
         for (i, transform) in self.forward_transforms.iter().enumerate() {
             for r in result.lanes_mut(Axis_nd(i + offset)).into_iter() {
@@ -211,7 +223,7 @@ impl<T: FftNum, D: Dimension> PeriodicConvolver<T, D> {
     }
 
     fn inverse_transform<D2: Dimension>(&self, mut f: Array<Complex<T>, D2>) -> Array<T, D2> {
-        let offset = D2::NDIM.unwrap() - D::NDIM.unwrap();
+        let offset = f.ndim() - self.k_abs.ndim();
         for (i, transform) in self.inverse_transforms.iter().enumerate() {
             for r in f.lanes_mut(Axis_nd(i + offset)).into_iter() {
                 self.transform(transform, r);

diff --git a/crates/feos-dft/src/functional.rs b/crates/feos-dft/src/functional.rs
@@ -11,6 +11,7 @@ use petgraph::graph::{Graph, UnGraph};
 use petgraph::visit::EdgeRef;
 use std::borrow::Cow;
 use std::ops::{Deref, MulAssign};
+use std::time::{Duration, Instant};
 
 impl<I: Clone, F: HelmholtzEnergyFunctionalDyn> HelmholtzEnergyFunctionalDyn
     for EquationOfState<Vec<I>, F>
@@ -112,17 +113,23 @@ 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>,
-    ) -> FeosResult<(Array<N, D>, Array<N, D::Larger>)>
+    ) -> FeosResult<(Array<N, D>, Array<N, D::Larger>, [Duration; 3])>
     where
         D: Dimension,
         D::Larger: Dimension<Smaller = D>,
     {
+        // calculate weighted densities
+        let start = Instant::now();
         let weighted_densities = convolver.weighted_densities(density);
+        let t_wd = start.elapsed();
+
+        // calculate partial derivatives
+        let start = Instant::now();
         let contributions = self.contributions();
         let mut partial_derivatives = Vec::new();
         let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
@@ -140,9 +147,17 @@ pub trait HelmholtzEnergyFunctional: Residual {
             partial_derivatives.push(pd);
             helmholtz_energy_density += &phi;
         }
+        let t_pd = start.elapsed();
+
+        // calculate functional derivative
+        let start = Instant::now();
+        let functional_derivative = convolver.functional_derivative(&partial_derivatives);
+        let t_fd = start.elapsed();
+
         Ok((
             helmholtz_energy_density,
-            convolver.functional_derivative(&partial_derivatives),
+            functional_derivative,
+            [t_wd, t_pd, t_fd],
         ))
     }
 

diff --git a/crates/feos-dft/src/functional_contribution.rs b/crates/feos-dft/src/functional_contribution.rs
@@ -1,10 +1,17 @@
 use crate::weight_functions::WeightFunctionInfo;
 use feos_core::{FeosResult, StateHD};
 use ndarray::RemoveAxis;
+#[cfg(feature = "rayon")]
+use ndarray::parallel::prelude::*;
 use ndarray::prelude::*;
 use num_dual::*;
 use num_traits::Zero;
 
+// The chunk size is specifically chosen to be the default amount of grid points in a 1D density
+// profile to switch to parallel evaluation only for 2D and 3D profiles.
+#[cfg(feature = "rayon")]
+const CHUNK_SIZE: usize = 2048;
+
 /// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers.
 pub trait FunctionalContribution: Sync + Send {
     /// Return the name of the contribution.
@@ -31,6 +38,36 @@ pub trait FunctionalContribution: Sync + Send {
         weighted_densities: ArrayView2<N>,
     ) -> FeosResult<Array1<N>>;
 
+    fn helmholtz_energy_density_parallel<N: DualNumCopy<Primitive = f64>>(
+        &self,
+        temperature: N,
+        weighted_densities: ArrayView2<N>,
+    ) -> FeosResult<Array1<N>> {
+        #[cfg(feature = "rayon")]
+        {
+            let n = weighted_densities.shape()[1];
+            if n > CHUNK_SIZE {
+                let mut phi = Array::zeros(n);
+                let mut phi_view = phi.view_mut();
+                let wd_iter = weighted_densities
+                    .axis_chunks_iter(Axis(1), CHUNK_SIZE)
+                    .into_par_iter();
+                let phi_iter = phi_view
+                    .axis_chunks_iter_mut(Axis(0), CHUNK_SIZE)
+                    .into_par_iter();
+                wd_iter.zip(phi_iter).try_for_each(|(wd, mut phi)| {
+                    self.helmholtz_energy_density_parallel(temperature, wd.view())
+                        .map(|p| phi.assign(&p))
+                })?;
+                Ok(phi)
+            } else {
+                self.helmholtz_energy_density(temperature, weighted_densities)
+            }
+        }
+        #[cfg(not(feature = "rayon"))]
+        self.helmholtz_energy_density(temperature, weighted_densities)
+    }
+
     fn bulk_helmholtz_energy_density<N: DualNum<Primitive = f64> + Copy>(
         &self,
         state: &StateHD<N>,
@@ -52,7 +89,7 @@ pub trait FunctionalContribution: Sync + Send {
             .unwrap()[0]
     }
 
-    fn first_partial_derivatives<N: DualNum<Primitive = f64> + Copy>(
+    fn first_partial_derivatives<N: DualNumCopy<Primitive = f64>>(
         &self,
         temperature: N,
         weighted_densities: Array2<N>,
@@ -66,7 +103,7 @@ pub trait FunctionalContribution: Sync + Send {
         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())?;
+            phi = self.helmholtz_energy_density_parallel(t, wd.view())?;
             first_partial_derivative
                 .index_axis_mut(Axis(0), i)
                 .assign(&phi.mapv(|p| p.eps));
@@ -93,7 +130,7 @@ pub trait FunctionalContribution: Sync + Send {
             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())?;
+                phi = self.helmholtz_energy_density_parallel(t, wd.view())?;
                 let p = phi.mapv(|p| p.eps1eps2);
                 second_partial_derivative
                     .index_axis_mut(Axis(0), i)

diff --git a/crates/feos-dft/src/geometry.rs b/crates/feos-dft/src/geometry.rs
@@ -6,7 +6,9 @@ use std::f64::consts::{FRAC_PI_3, PI};
 /// Grids with up to three dimensions.
 #[derive(Clone)]
 pub enum Grid {
+    Bulk,
     Cartesian1(Axis),
+    Periodical1(Axis),
     Cartesian2(Axis, Axis),
     Periodical2(Axis, Axis, Angle),
     Cartesian3(Axis, Axis, Axis),
@@ -27,7 +29,8 @@ impl Grid {
 
     pub fn axes(&self) -> Vec<&Axis> {
         match self {
-            Self::Cartesian1(x) => vec![x],
+            Self::Bulk => vec![],
+            Self::Cartesian1(x) | Self::Periodical1(x) => vec![x],
             Self::Cartesian2(x, y) | Self::Periodical2(x, y, _) => vec![x, y],
             Self::Cartesian3(x, y, z) | Self::Periodical3(x, y, z, _) => vec![x, y, z],
             Self::Spherical(r) | Self::Polar(r) => vec![r],
@@ -37,7 +40,8 @@ impl Grid {
 
     pub fn axes_mut(&mut self) -> Vec<&mut Axis> {
         match self {
-            Self::Cartesian1(x) => vec![x],
+            Self::Bulk => vec![],
+            Self::Cartesian1(x) | Self::Periodical1(x) => vec![x],
             Self::Cartesian2(x, y) | Self::Periodical2(x, y, _) => vec![x, y],
             Self::Cartesian3(x, y, z) | Self::Periodical3(x, y, z, _) => vec![x, y, z],
             Self::Spherical(r) | Self::Polar(r) => vec![r],
@@ -72,7 +76,11 @@ impl Grid {
 
     pub fn mesh(&self) -> Vec<Length<ArrayD<f64>>> {
         match self {
-            Grid::Cartesian1(ax) | Grid::Spherical(ax) | Grid::Polar(ax) => {
+            Self::Bulk => vec![],
+            Grid::Cartesian1(ax)
+            | Self::Periodical1(ax)
+            | Grid::Spherical(ax)
+            | Grid::Polar(ax) => {
                 vec![Length::from_reduced(ax.grid.clone()).into_dyn()]
             }
             Grid::Cartesian2(u, v) => mesh_2d(u, v, 90.0 * DEGREES),