prehner · Oct 26, 2026 · Oct 16, 2026 · Oct 16, 2026 · Oct 16, 2026 · Oct 21, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Changed
+- Removed most heap allocations from the evaluation of the residual Helmholtz energy of PC-SAFT. [#397](https://github.com/feos-org/feos/pull/397)
+- `ResidualDyn::for_each_contribution` and `Residual::for_each_contribution`, which pass each Helmholtz energy contribution to a closure, are now the only methods a model has to implement to evaluate the Helmholtz energy. `Residual::reduced_residual_helmholtz_energy_density` and `Residual::reduced_helmholtz_energy_density_contributions` are provided methods built on top of it, and both methods were removed from `ResidualDyn`. [#397](https://github.com/feos-org/feos/pull/397)
+- `HelmholtzEnergyFunctional::evaluate_bulk` passes each contribution to a closure instead of returning a `Vec`. [#397](https://github.com/feos-org/feos/pull/397)
+- `MonomerShape::NonSpherical` holds a reference `&DVector<f64>` to the segment numbers instead of an owned `DVector<D>`. [#397](https://github.com/feos-org/feos/pull/397)
+- `StateHD::new_density` takes the partial densities by value. [#397](https://github.com/feos-org/feos/pull/397)
 
 ## [0.10.4] - 2026-10-21
 ### Fixed

diff --git a/crates/feos-core/src/cubic.rs b/crates/feos-core/src/cubic.rs
@@ -125,10 +125,11 @@ impl ResidualDyn for PengRobinson {
         D::from(0.9) / molefracs.dot(&self.b.map(D::from))
     }
 
-    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
+    fn for_each_contribution<D: DualNum<f64> + Copy>(
         &self,
         state: &StateHD<D>,
-    ) -> Vec<(&'static str, D)> {
+        mut f: impl FnMut(&'static str, D),
+    ) {
         let density = state.partial_density.sum();
         let x = &state.molefracs;
         let ak = &self
@@ -149,11 +150,11 @@ impl ResidualDyn for PengRobinson {
 
         // Helmholtz energy
         let v = density.recip();
-        let f = density
+        let a = density
             * ((v / (v - b)).ln()
                 - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
                     * ((v + b * (1.0 + SQRT_2)) / (v + b * (1.0 - SQRT_2))).ln());
-        vec![("Peng Robinson", f)]
+        f("Peng Robinson", a)
     }
 }
 

diff --git a/crates/feos-core/src/equation_of_state/mod.rs b/crates/feos-core/src/equation_of_state/mod.rs
@@ -57,12 +57,12 @@ impl<I, R: ResidualDyn> ResidualDyn for EquationOfState<Vec<I>, R> {
         self.residual.compute_max_density(molefracs)
     }
 
-    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
+    fn for_each_contribution<D: DualNum<f64> + Copy>(
         &self,
         state: &StateHD<D>,
-    ) -> Vec<(&'static str, D)> {
-        self.residual
-            .reduced_helmholtz_energy_density_contributions(state)
+        f: impl FnMut(&'static str, D),
+    ) {
+        self.residual.for_each_contribution(state, f)
     }
 }
 
@@ -99,17 +99,8 @@ impl<I: Clone, R: Residual<Const<N>, D>, D: DualNum<f64> + Copy, const N: usize>
         self.residual.compute_max_density(molefracs)
     }
 
-    fn reduced_helmholtz_energy_density_contributions(
-        &self,
-        state: &StateHD<D, Const<N>>,
-    ) -> Vec<(&'static str, D)> {
-        self.residual
-            .reduced_helmholtz_energy_density_contributions(state)
-    }
-
-    fn reduced_residual_helmholtz_energy_density(&self, state: &StateHD<D, Const<N>>) -> D {
-        self.residual
-            .reduced_residual_helmholtz_energy_density(state)
+    fn for_each_contribution(&self, state: &StateHD<D, Const<N>>, f: impl FnMut(&'static str, D)) {
+        self.residual.for_each_contribution(state, f)
     }
 }
 

diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs
@@ -65,11 +65,12 @@ pub trait ResidualDyn {
     fn compute_max_density<D: DualNum<f64> + Copy>(&self, molefracs: &DVector<D>) -> D;
 
     /// Evaluate the reduced Helmholtz energy density of each individual contribution
-    /// and return them together with a string representation of the contribution.
-    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
+    /// and pass it to `f` together with a string representation of the contribution.
+    fn for_each_contribution<D: DualNum<f64> + Copy>(
         &self,
         state: &StateHD<D>,
-    ) -> Vec<(&'static str, D)>;
+        f: impl FnMut(&'static str, D),
+    );
 }
 
 impl<C: Deref<Target = T> + Clone, T: ResidualDyn, D: DualNum<f64> + Copy> Residual<Dyn, D> for C {
@@ -87,11 +88,8 @@ impl<C: Deref<Target = T> + Clone, T: ResidualDyn, D: DualNum<f64> + Copy> Resid
     fn compute_max_density(&self, molefracs: &DVector<D>) -> D {
         ResidualDyn::compute_max_density(self.deref(), molefracs)
     }
-    fn reduced_helmholtz_energy_density_contributions(
-        &self,
-        state: &StateHD<D, Dyn>,
-    ) -> Vec<(&'static str, D)> {
-        ResidualDyn::reduced_helmholtz_energy_density_contributions(self.deref(), state)
+    fn for_each_contribution(&self, state: &StateHD<D, Dyn>, f: impl FnMut(&'static str, D)) {
+        ResidualDyn::for_each_contribution(self.deref(), state, f)
     }
 }
 
@@ -130,18 +128,26 @@ where
     /// equation of state anyways).
     fn compute_max_density(&self, molefracs: &OVector<D, N>) -> D;
 
+    /// Evaluate the reduced Helmholtz energy density of each individual contribution
+    /// and pass it to `f` together with a string representation of the contribution.
+    fn for_each_contribution(&self, state: &StateHD<D, N>, f: impl FnMut(&'static str, D));
+
     /// Evaluate the reduced Helmholtz energy density of each individual contribution
     /// and return them together with a string representation of the contribution.
     fn reduced_helmholtz_energy_density_contributions(
         &self,
         state: &StateHD<D, N>,
-    ) -> Vec<(&'static str, D)>;
+    ) -> Vec<(&'static str, D)> {
+        let mut v = Vec::new();
+        self.for_each_contribution(state, |name, a| v.push((name, a)));
+        v
+    }
 
     /// Evaluate the residual reduced Helmholtz energy density $\beta f^\mathrm{res}$.
     fn reduced_residual_helmholtz_energy_density(&self, state: &StateHD<D, N>) -> D {
-        self.reduced_helmholtz_energy_density_contributions(state)
-            .iter()
-            .fold(D::zero(), |acc, (_, a)| acc + a)
+        let mut a = D::zero();
+        self.for_each_contribution(state, |_, a_i| a += a_i);
+        a
     }
 
     /// Evaluate the Helmholtz energy of each individual contribution and return them
@@ -153,10 +159,9 @@ where
         moles: &OVector<D, N>,
     ) -> Vec<(&'static str, D)> {
         let state = StateHD::new(temperature, volume, moles);
-        self.reduced_helmholtz_energy_density_contributions(&state)
-            .into_iter()
-            .map(|(n, f)| (n, f * temperature * volume))
-            .collect()
+        let mut v = Vec::new();
+        self.for_each_contribution(&state, |name, a| v.push((name, a * temperature * volume)));
+        v
     }
 
     /// Evaluate the residual Helmholtz energy $A^\mathrm{res}$.
@@ -360,7 +365,7 @@ where
     {
         let (f_res, mu_res, dmu_res) = N::hessian(
             |rho, &t| {
-                let state = StateHD::new_density(t, &rho);
+                let state = StateHD::new_density(t, rho);
                 self.lift()
                     .reduced_residual_helmholtz_energy_density(&state)
                     * t
@@ -403,7 +408,7 @@ where
     {
         let (_, _, f_res_t, mu_res_t) = N::partial_hessian(
             |rho, t, _: &()| {
-                let state = StateHD::new_density(t, &rho);
+                let state = StateHD::new_density(t, rho);
                 self.lift()
                     .reduced_residual_helmholtz_energy_density(&state)
                     * t
@@ -506,10 +511,10 @@ impl ResidualDyn for NoResidual {
         D::one()
     }
 
-    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
+    fn for_each_contribution<D: DualNum<f64> + Copy>(
         &self,
         _: &StateHD<D>,
-    ) -> Vec<(&'static str, D)> {
-        vec![]
+        _: impl FnMut(&'static str, D),
+    ) {
     }
 }
diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs
@@ -949,6 +949,9 @@ fn log_iteration(
     x2: &[f64],
     newton: bool,
 ) {
+    if verbosity < Verbosity::Iter {
+        return;
+    }
     let error = error.map_or_else(|| format!("{:14}", ""), |e| format!("{:<14.8e}", e.abs()));
     log_iter!(
         verbosity,

diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs
@@ -82,17 +82,16 @@ where
 {
     /// Create a new `StateHD` for given temperature, volume and composition.
     pub fn new(temperature: D, volume: D, moles: &OVector<D, N>) -> Self {
-        Self::new_density(temperature, &(moles / volume))
+        Self::new_density(temperature, moles / volume)
     }
 
     /// Create a new `StateHD` for given temperature and partial densities
-    pub fn new_density(temperature: D, partial_density: &OVector<D, N>) -> Self {
-        let molefracs = partial_density / partial_density.sum();
-
+    pub fn new_density(temperature: D, partial_density: OVector<D, N>) -> Self {
+        let molefracs = &partial_density / partial_density.sum();
         Self {
             temperature,
             molefracs,
-            partial_density: partial_density.clone(),
+            partial_density,
         }
     }
 

diff --git a/crates/feos-derive/src/residual.rs b/crates/feos-derive/src/residual.rs
@@ -36,10 +36,10 @@ fn impl_residual(
             Self::#name(residual) => residual.compute_max_density(moles)
         }
     });
-    let reduced_helmholtz_energy_density_contributions = variants.iter().map(|v| {
+    let for_each_contribution = variants.iter().map(|v| {
         let name = &v.ident;
         quote! {
-            Self::#name(residual) => residual.reduced_helmholtz_energy_density_contributions(state)
+            Self::#name(residual) => residual.for_each_contribution(state, f)
         }
     });
     // let subset = variants.iter().map(|v| {
@@ -61,9 +61,9 @@ fn impl_residual(
                     #(#compute_max_density,)*
                 }
             }
-            fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> Vec<(&'static str, D)> {
+            fn for_each_contribution<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>, f: impl FnMut(&'static str, D)) {
                 match self {
-                    #(#reduced_helmholtz_energy_density_contributions,)*
+                    #(#for_each_contribution,)*
                 }
             }
             // fn subset(&self, component_list: &[usize]) -> Self {

diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs
@@ -301,11 +301,12 @@ impl ResidualDyn for Helium {
         D::from(1.0)
     }
 
-    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
+    fn for_each_contribution<D: DualNum<f64> + Copy>(
         &self,
         state: &StateHD<D>,
-    ) -> Vec<(&'static str, D)> {
-        self.evaluate_bulk(state)
+        f: impl FnMut(&'static str, D),
+    ) {
+        self.evaluate_bulk(state, f)
     }
 }
 

diff --git a/crates/feos-dft/src/functional.rs b/crates/feos-dft/src/functional.rs
@@ -228,17 +228,19 @@ pub trait HelmholtzEnergyFunctional: Residual {
         i
     }
 
-    fn evaluate_bulk<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> Vec<(&'static str, D)> {
-        let mut res: Vec<_> = self
-            .contributions()
-            .map(|c| (c.name(), c.bulk_helmholtz_energy_density(state)))
-            .collect();
-        res.push((
+    fn evaluate_bulk<D: DualNum<f64> + Copy>(
+        &self,
+        state: &StateHD<D>,
+        mut f: impl FnMut(&'static str, D),
+    ) {
+        for c in self.contributions() {
+            f(c.name(), c.bulk_helmholtz_energy_density(state))
+        }
+        f(
             "Ideal chain",
             self.ideal_chain_contribution()
                 .bulk_helmholtz_energy_density(&state.partial_density),
-        ));
-        res
+        )
     }
 }
 

diff --git a/crates/feos/benches/contributions.rs b/crates/feos/benches/contributions.rs
@@ -74,13 +74,14 @@ fn pcsaft(c: &mut Criterion) {
                 State::new_npt(&&eos, t, p, &moles, Some(DensityInitialization::Liquid)).unwrap();
             let temperature = Dual64::from(state.temperature.into_reduced()).derivative();
             let molar_volume = Dual::from(1.0 / state.density.into_reduced());
-            let moles = state.moles().unwrap().to_reduced().map(Dual::from);
+            // one particle in the volume per particle, i.e., at the state's density
+            let molefracs = state.molefracs.map(Dual::from);
             // let state_hd = state.derive1(Derivative::DT);
             let name1 = comp1.identifier.name.as_deref().unwrap();
             let name2 = comp2.identifier.name.as_deref().unwrap();
             let mix = format!("{name1}_{name2}");
             group.bench_function(mix, |b| {
-                b.iter(|| (&eos).residual_helmholtz_energy(temperature, molar_volume, &moles))
+                b.iter(|| (&eos).residual_helmholtz_energy(temperature, molar_volume, &molefracs))
             });
         }
     }

diff --git a/crates/feos/benches/dual_numbers_saftvrmie.rs b/crates/feos/benches/dual_numbers_saftvrmie.rs
@@ -23,7 +23,7 @@ fn state_saftvrmie(n: usize, eos: &SaftVRMie) -> State<&SaftVRMie> {
     State::new_nvt(
         &eos,
         temperature,
-        cp.volume().unwrap(),
+        10.0 * MOL / cp.density,
         &(molefracs * 10. * MOL),
     )
     .unwrap()

diff --git a/crates/feos/benches/state_creation.rs b/crates/feos/benches/state_creation.rs
@@ -113,9 +113,6 @@ fn bench_states<E: Residual>(c: &mut Criterion, group_name: &str, eos: &E) {
         group.bench_function("critical_point_binary_t", |b| {
             b.iter(|| critical_point_binary((eos, crit.temperature)))
         });
-        group.bench_function("critical_point_binary_p", |b| {
-            b.iter(|| critical_point_binary((eos, crit.pressure(Contributions::Total))))
-        });
     }
     if ncomponents != 1 {
         group.bench_function("tp_flash", |b| {