finalized adaption of PC-SAFT family

feos-org · g-bauer · Oct 13, 2025 · Aug 18, 2025 · Aug 18, 2025 · Aug 18, 2025
commit 2826fa434f3d90c1985aef879770980a966cab07
diff --git a/crates/feos-core/src/state/residual_properties.rs b/crates/feos-core/src/state/residual_properties.rs
@@ -578,7 +578,7 @@ where
     /// Pressure $p$ evaluated for each contribution of the equation of state.
     pub fn pressure_contributions(&self) -> Vec<(String, Pressure<D>)> {
         let t = Dual::from_re(self.temperature.into_reduced());
-        let v = Dual::from_re(self.temperature.into_reduced()).derivative();
+        let v = Dual::from_re(self.density.into_reduced().recip()).derivative();
         let x = self.molefracs.map(Dual::from_re);
         let contributions = self
             .eos

diff --git a/crates/feos-dft/src/functional.rs b/crates/feos-dft/src/functional.rs
@@ -105,21 +105,23 @@ pub trait HelmholtzEnergyFunctional: Residual {
             .collect()
     }
 
-    fn m(&self) -> Cow<'_, DVector<f64>> {
+    fn m(&self) -> Cow<'_, [f64]> {
         match self.molecule_shape() {
-            MoleculeShape::Spherical(n) => Cow::Owned(DVector::from_element(n, 1.0)),
-            MoleculeShape::NonSpherical(m) => Cow::Borrowed(m),
+            MoleculeShape::Spherical(n) => Cow::Owned(vec![1.0; n]),
+            MoleculeShape::NonSpherical(m) => Cow::Borrowed(m.as_slice()),
             MoleculeShape::Heterosegmented(component_index) => {
-                Cow::Owned(DVector::from_element(component_index.len(), 1.0))
+                Cow::Owned(vec![1.0; component_index.len()])
             }
         }
     }
 
-    fn component_index(&self) -> Cow<'_, DVector<usize>> {
+    fn component_index(&self) -> Cow<'_, [usize]> {
         match self.molecule_shape() {
-            MoleculeShape::Spherical(n) => Cow::Owned(DVector::from_fn(n, |i, _| i)),
-            MoleculeShape::NonSpherical(m) => Cow::Owned(DVector::from_fn(m.len(), |i, _| i)),
-            MoleculeShape::Heterosegmented(component_index) => Cow::Borrowed(component_index),
+            MoleculeShape::Spherical(n) => Cow::Owned((0..n).collect()),
+            MoleculeShape::NonSpherical(m) => Cow::Owned((0..m.len()).collect()),
+            MoleculeShape::Heterosegmented(component_index) => {
+                Cow::Borrowed(component_index.as_slice())
+            }
         }
     }
 

diff --git a/crates/feos-dft/src/ideal_chain_contribution.rs b/crates/feos-dft/src/ideal_chain_contribution.rs
@@ -11,10 +11,10 @@ pub struct IdealChainContribution {
 }
 
 impl IdealChainContribution {
-    pub fn new(component_index: &DVector<usize>, m: &DVector<f64>) -> Self {
+    pub fn new(component_index: &[usize], m: &[f64]) -> Self {
         Self {
-            component_index: component_index.clone(),
-            m: m.clone(),
+            component_index: DVector::from_column_slice(component_index),
+            m: DVector::from_column_slice(m),
         }
     }
 

diff --git a/crates/feos-dft/src/weight_functions.rs b/crates/feos-dft/src/weight_functions.rs
@@ -169,7 +169,7 @@ pub enum WeightFunctionShape {
 /// Information about weight functions
 pub struct WeightFunctionInfo<T> {
     /// Index of the component that each individual segment belongs to.
-    pub(crate) component_index: Vec<usize>,
+    pub(crate) component_index: DVector<usize>,
     /// Flag if local density is required in the functional
     pub(crate) local_density: bool,
     /// Container for scalar component-wise weighted densities
@@ -197,7 +197,7 @@ impl<T> WeightFunctionInfo<T> {
 
 impl<T> WeightFunctionInfo<T> {
     /// Initializing empty `WeightFunctionInfo`.
-    pub fn new(component_index: Vec<usize>, local_density: bool) -> Self {
+    pub fn new(component_index: DVector<usize>, local_density: bool) -> Self {
         Self {
             component_index,
             local_density,

diff --git a/crates/feos/src/association/dft.rs b/crates/feos/src/association/dft.rs
@@ -42,7 +42,7 @@ where
         let p = self.model;
         let r = p.hs_diameter(temperature) * N::from(0.5);
         let [_, _, _, c3] = p.geometry_coefficients(temperature);
-        WeightFunctionInfo::new(p.component_index().into_owned(), false)
+        WeightFunctionInfo::new(p.component_index().into_owned().into(), false)
             .add(
                 WeightFunction::new_scaled(r.clone(), WeightFunctionShape::Delta),
                 false,

diff --git a/crates/feos/src/epcsaft/eos/born.rs b/crates/feos/src/epcsaft/eos/born.rs
@@ -2,17 +2,17 @@ use super::ElectrolytePcSaftVariants;
 use crate::epcsaft::eos::permittivity::Permittivity;
 use crate::epcsaft::parameters::ElectrolytePcSaftPars;
 use feos_core::StateHD;
-use ndarray::Array1;
+use nalgebra::DVector;
 use num_dual::DualNum;
 
 pub struct Born;
 
 impl Born {
-    pub fn helmholtz_energy<D: DualNum<f64> + Copy>(
+    pub fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
         &self,
         parameters: &ElectrolytePcSaftPars,
         state: &StateHD<D>,
-        diameter: &Array1<D>,
+        diameter: &DVector<D>,
     ) -> D {
         // Parameters
         let p = parameters;
@@ -25,14 +25,14 @@ impl Born {
             .unwrap()
             .permittivity;
 
-        // Calculate sum xi zi^2 / di
-        let mut sum_xi_zi_ai = D::zero();
+        // Calculate sum rhoi zi^2 / di
+        let mut sum_rhoi_zi_ai = D::zero();
         for i in 0..state.molefracs.len() {
-            sum_xi_zi_ai += state.molefracs[i] * p.z[i].powi(2) / diameter[i];
+            sum_rhoi_zi_ai += state.partial_density[i] * p.z[i].powi(2) / diameter[i];
         }
 
         // Calculate born contribution
-        -lambda_b * (epsilon_r - 1.) * sum_xi_zi_ai * state.moles.sum()
+        -lambda_b * (epsilon_r - 1.) * sum_rhoi_zi_ai
     }
 }
 
@@ -42,17 +42,17 @@ mod tests {
     use crate::epcsaft::parameters::utils::{water_nacl_parameters, water_nacl_parameters_perturb};
     use crate::hard_sphere::HardSphereProperties;
     use approx::assert_relative_eq;
-    use ndarray::arr1;
+    use nalgebra::dvector;
 
     #[test]
     fn helmholtz_energy_perturb() {
         let p = ElectrolytePcSaftPars::new(&water_nacl_parameters_perturb()).unwrap();
 
         let t = 298.0;
         let v = 31.875;
-        let s = StateHD::new(t, v, arr1(&[0.9, 0.05, 0.05]));
+        let s = StateHD::new(t, v, &dvector![0.9, 0.05, 0.05]);
         let d = p.hs_diameter(t);
-        let a_rust = Born.helmholtz_energy(&p, &s, &d);
+        let a_rust = Born.helmholtz_energy_density(&p, &s, &d) * v;
 
         assert_relative_eq!(a_rust, -22.51064553710294, epsilon = 1e-10);
     }
@@ -63,9 +63,9 @@ mod tests {
 
         let t = 298.0;
         let v = 31.875;
-        let s = StateHD::new(t, v, arr1(&[0.9, 0.05, 0.05]));
+        let s = StateHD::new(t, v, &dvector![0.9, 0.05, 0.05]);
         let d = p.hs_diameter(t);
-        let a_rust = Born.helmholtz_energy(&p, &s, &d);
+        let a_rust = Born.helmholtz_energy_density(&p, &s, &d) * v;
 
         assert_relative_eq!(a_rust, -22.525624511559244, epsilon = 1e-10);
     }

diff --git a/crates/feos/src/epcsaft/eos/dispersion.rs b/crates/feos/src/epcsaft/eos/dispersion.rs
@@ -1,8 +1,8 @@
 use crate::epcsaft::parameters::ElectrolytePcSaftPars;
 use feos_core::StateHD;
-use ndarray::{Array, Array1, Array2};
+use nalgebra::{DMatrix, DVector};
 use num_dual::DualNum;
-use std::f64::consts::{FRAC_PI_3, PI};
+use std::f64::consts::{FRAC_PI_6, PI};
 
 pub const A0: [f64; 7] = [
     0.91056314451539,
@@ -62,34 +62,34 @@ pub const B2: [f64; 7] = [
 pub const T_REF: f64 = 298.15;
 
 impl ElectrolytePcSaftPars {
-    pub fn k_ij_t<D: DualNum<f64>>(&self, temperature: D) -> Array2<f64> {
+    pub fn k_ij_t<D: DualNum<f64>>(&self, temperature: D) -> DMatrix<f64> {
         let k_ij = &self.k_ij;
         let n = self.m.len();
 
-        let mut k_ij_t = Array::zeros((n, n));
+        let mut k_ij_t = DMatrix::zeros(n, n);
 
         for i in 0..n {
             for j in 0..n {
                 // Calculate k_ij(T)
-                k_ij_t[[i, j]] = (temperature.re() - T_REF) * k_ij[[i, j]][1]
-                    + (temperature.re() - T_REF).powi(2) * k_ij[[i, j]][2]
-                    + (temperature.re() - T_REF).powi(3) * k_ij[[i, j]][3]
-                    + k_ij[[i, j]][0];
+                k_ij_t[(i, j)] = (temperature.re() - T_REF) * k_ij[(i, j)][1]
+                    + (temperature.re() - T_REF).powi(2) * k_ij[(i, j)][2]
+                    + (temperature.re() - T_REF).powi(3) * k_ij[(i, j)][3]
+                    + k_ij[(i, j)][0];
             }
         }
         //println!("k_ij_t: {}", k_ij_t);
         k_ij_t
     }
 
-    pub fn epsilon_k_ij_t<D: DualNum<f64>>(&self, temperature: D) -> Array2<f64> {
+    pub fn epsilon_k_ij_t<D: DualNum<f64>>(&self, temperature: D) -> DMatrix<f64> {
         let k_ij_t = self.k_ij_t(temperature);
         let n = self.m.len();
 
-        let mut epsilon_k_ij_t = Array::zeros((n, n));
+        let mut epsilon_k_ij_t = DMatrix::zeros(n, n);
 
         for i in 0..n {
             for j in 0..n {
-                epsilon_k_ij_t[[i, j]] = (1.0 - k_ij_t[[i, j]]) * self.e_k_ij[[i, j]];
+                epsilon_k_ij_t[(i, j)] = (1.0 - k_ij_t[(i, j)]) * self.e_k_ij[(i, j)];
             }
         }
         epsilon_k_ij_t
@@ -99,31 +99,28 @@ impl ElectrolytePcSaftPars {
 pub struct Dispersion;
 
 impl Dispersion {
-    pub fn helmholtz_energy<D: DualNum<f64> + Copy>(
+    pub fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
         &self,
         parameters: &ElectrolytePcSaftPars,
         state: &StateHD<D>,
-        diameter: &Array1<D>,
+        diameter: &DVector<D>,
     ) -> D {
         // auxiliary variables
         let n = parameters.m.len();
         let p = parameters;
         let rho = &state.partial_density;
 
-        // temperature dependent segment radius
-        let r = diameter * 0.5;
-
         // convert sigma_ij
         let sigma_ij_t = p.sigma_ij_t(state.temperature);
 
         // convert epsilon_k_ij
         let epsilon_k_ij_t = p.epsilon_k_ij_t(state.temperature);
 
         // packing fraction
-        let eta = (rho * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3;
+        let eta = rho.dot(&diameter.zip_map(&p.m, |d, m| d.powi(3) * m)) * FRAC_PI_6;
 
         // mean segment number
-        let m = (&state.molefracs * &p.m).sum();
+        let m = state.molefracs.dot(&p.m.map(D::from));
 
         // mixture densities, crosswise interactions of all segments on all chains
         let mut rho1mix = D::zero();
@@ -154,7 +151,7 @@ impl Dispersion {
             .recip();
 
         // Helmholtz energy
-        (-rho1mix * i1 * 2.0 - rho2mix * m * c1 * i2) * PI * state.volume
+        (-rho1mix * i1 * 2.0 - rho2mix * m * c1 * i2) * PI
     }
 }
 
@@ -166,17 +163,17 @@ mod tests {
     };
     use crate::hard_sphere::HardSphereProperties;
     use approx::assert_relative_eq;
-    use ndarray::arr1;
+    use nalgebra::dvector;
 
     #[test]
     fn helmholtz_energy() {
         let p = ElectrolytePcSaftPars::new(&propane_parameters()).unwrap();
         let t = 250.0;
         let v = 1000.0;
         let n = 1.0;
-        let s = StateHD::new(t, v, arr1(&[n]));
+        let s = StateHD::new(t, v, &dvector![n]);
         let d = p.hs_diameter(t);
-        let a_rust = Dispersion.helmholtz_energy(&p, &s, &d);
+        let a_rust = Dispersion.helmholtz_energy_density(&p, &s, &d) * v;
         assert_relative_eq!(a_rust, -1.0622531100351962, epsilon = 1e-10);
     }
 
@@ -188,13 +185,13 @@ mod tests {
         let t = 250.0;
         let v = 2.5e28;
         let n = 1.0;
-        let s = StateHD::new(t, v, arr1(&[n]));
-        let a1 = Dispersion.helmholtz_energy(&p1, &s, &p1.hs_diameter(t));
-        let a2 = Dispersion.helmholtz_energy(&p2, &s, &p2.hs_diameter(t));
-        let s1m = StateHD::new(t, v, arr1(&[n, 0.0]));
-        let a1m = Dispersion.helmholtz_energy(&p12, &s1m, &p12.hs_diameter(t));
-        let s2m = StateHD::new(t, v, arr1(&[0.0, n]));
-        let a2m = Dispersion.helmholtz_energy(&p12, &s2m, &p12.hs_diameter(t));
+        let s = StateHD::new(t, v, &dvector![n]);
+        let a1 = Dispersion.helmholtz_energy_density(&p1, &s, &p1.hs_diameter(t));
+        let a2 = Dispersion.helmholtz_energy_density(&p2, &s, &p2.hs_diameter(t));
+        let s1m = StateHD::new(t, v, &dvector![n, 0.0]);
+        let a1m = Dispersion.helmholtz_energy_density(&p12, &s1m, &p12.hs_diameter(t));
+        let s2m = StateHD::new(t, v, &dvector![0.0, n]);
+        let a2m = Dispersion.helmholtz_energy_density(&p12, &s2m, &p12.hs_diameter(t));
         assert_relative_eq!(a1, a1m, epsilon = 1e-14);
         assert_relative_eq!(a2, a2m, epsilon = 1e-14);
     }

diff --git a/crates/feos/src/epcsaft/eos/hard_chain.rs b/crates/feos/src/epcsaft/eos/hard_chain.rs
@@ -1,14 +1,14 @@
 use crate::epcsaft::parameters::ElectrolytePcSaftPars;
 use crate::hard_sphere::HardSphereProperties;
 use feos_core::StateHD;
-use ndarray::Array;
+use nalgebra::DVector;
 use num_dual::*;
 
 pub struct HardChain;
 
 impl HardChain {
     #[inline]
-    pub fn helmholtz_energy<D: DualNum<f64> + Copy>(
+    pub fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
         &self,
         parameters: &ElectrolytePcSaftPars,
         state: &StateHD<D>,
@@ -18,13 +18,11 @@ impl HardChain {
         let [zeta2, zeta3] = p.zeta(state.temperature, &state.partial_density, [2, 3]);
         let frac_1mz3 = -(zeta3 - 1.0).recip();
         let c = zeta2 * frac_1mz3 * frac_1mz3;
-        let g_hs =
-            d.mapv(|d| frac_1mz3 + d * c * 1.5 - d.powi(2) * c.powi(2) * (zeta3 - 1.0) * 0.5);
-        Array::from_shape_fn(p.m.len(), |i| {
+        let g_hs = d.map(|d| frac_1mz3 + d * c * 1.5 - d.powi(2) * c.powi(2) * (zeta3 - 1.0) * 0.5);
+        DVector::from_fn(p.m.len(), |i, _| {
             state.partial_density[i] * (1.0 - p.m[i]) * g_hs[i].ln()
         })
         .sum()
-            * state.volume
     }
 }
 
@@ -35,16 +33,16 @@ mod tests {
         butane_parameters, propane_butane_parameters, propane_parameters,
     };
     use approx::assert_relative_eq;
-    use ndarray::arr1;
+    use nalgebra::dvector;
 
     #[test]
     fn helmholtz_energy() {
         let p = ElectrolytePcSaftPars::new(&propane_parameters()).unwrap();
         let t = 250.0;
         let v = 1000.0;
         let n = 1.0;
-        let s = StateHD::new(t, v, arr1(&[n]));
-        let a_rust = HardChain.helmholtz_energy(&p, &s);
+        let s = StateHD::new(t, v, &dvector![n]);
+        let a_rust = HardChain.helmholtz_energy_density(&p, &s) * v;
         assert_relative_eq!(a_rust, -0.12402626171926148, epsilon = 1e-10);
     }
 
@@ -56,13 +54,13 @@ mod tests {
         let t = 250.0;
         let v = 2.5e28;
         let n = 1.0;
-        let s = StateHD::new(t, v, arr1(&[n]));
-        let a1 = HardChain.helmholtz_energy(&p1, &s);
-        let a2 = HardChain.helmholtz_energy(&p2, &s);
-        let s1m = StateHD::new(t, v, arr1(&[n, 0.0]));
-        let a1m = HardChain.helmholtz_energy(&p12, &s1m);
-        let s2m = StateHD::new(t, v, arr1(&[0.0, n]));
-        let a2m = HardChain.helmholtz_energy(&p12, &s2m);
+        let s = StateHD::new(t, v, &dvector![n]);
+        let a1 = HardChain.helmholtz_energy_density(&p1, &s);
+        let a2 = HardChain.helmholtz_energy_density(&p2, &s);
+        let s1m = StateHD::new(t, v, &dvector![n, 0.0]);
+        let a1m = HardChain.helmholtz_energy_density(&p12, &s1m);
+        let s2m = StateHD::new(t, v, &dvector![0.0, n]);
+        let a2m = HardChain.helmholtz_energy_density(&p12, &s2m);
         assert_relative_eq!(a1, a1m, epsilon = 1e-14);
         assert_relative_eq!(a2, a2m, epsilon = 1e-14);
     }