feos-org · prehner · Jul 23, 2026
diff --git a/crates/feos-dft/src/adsorption/mod.rs b/crates/feos-dft/src/adsorption/mod.rs
@@ -17,7 +17,7 @@ mod fea_potential;
 mod pore;
 mod pore2d;
 pub use external_potential::{ExternalPotential, FluidParameters};
-pub use pore::{HenryCoefficient, Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
+pub use pore::{HenryCoefficient, Pore, Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
 pub use pore2d::{Pore2D, PoreProfile2D};
 
 #[cfg(feature = "rayon")]
@@ -54,7 +54,7 @@ where
     }
 
     /// Calculate an adsorption isotherm (starting at low pressure)
-    pub fn adsorption_isotherm<S: PoreSpecification<D>, X: Composition<f64, Dyn> + Clone>(
+    pub fn adsorption_isotherm<S: Pore<D>, X: Composition<f64, Dyn> + Clone>(
         functional: &F,
         temperature: Temperature,
         pressure: &Pressure<Array1<f64>>,
@@ -74,7 +74,7 @@ where
     }
 
     /// Calculate an desorption isotherm (starting at high pressure)
-    pub fn desorption_isotherm<S: PoreSpecification<D>, X: Composition<f64, Dyn> + Clone>(
+    pub fn desorption_isotherm<S: Pore<D>, X: Composition<f64, Dyn> + Clone>(
         functional: &F,
         temperature: Temperature,
         pressure: &Pressure<Array1<f64>>,
@@ -99,7 +99,7 @@ where
     }
 
     /// Calculate an equilibrium isotherm
-    pub fn equilibrium_isotherm<S: PoreSpecification<D>, X: Composition<f64, Dyn> + Clone>(
+    pub fn equilibrium_isotherm<S: Pore<D>, X: Composition<f64, Dyn> + Clone>(
         functional: &F,
         temperature: Temperature,
         pressure: &Pressure<Array1<f64>>,
@@ -182,7 +182,7 @@ where
         }
     }
 
-    fn isotherm<S: PoreSpecification<D>, X: Composition<f64, Dyn> + Clone>(
+    fn isotherm<S: Pore<D>, X: Composition<f64, Dyn> + Clone>(
         functional: &F,
         temperature: Temperature,
         pressure: &Pressure<Array1<f64>>,
@@ -210,7 +210,10 @@ where
                 _ => unreachable!(),
             };
         }
-        let profile = pore.initialize(&bulk, None, None)?.solve(solver)?.profile;
+        let profile = pore
+            .initialize(&bulk, None, None, PoreSpecification::ChemicalPotential)?
+            .solve(solver)?
+            .profile;
         let external_potential = Some(&profile.external_potential);
         let mut old_density = Some(&profile.density);
 
@@ -229,8 +232,18 @@ where
                     .clone();
             }
 
-            let p = pore.initialize(&bulk, old_density, external_potential)?;
-            let p2 = pore.initialize(&bulk, None, external_potential)?;
+            let p = pore.initialize(
+                &bulk,
+                old_density,
+                external_potential,
+                PoreSpecification::ChemicalPotential,
+            )?;
+            let p2 = pore.initialize(
+                &bulk,
+                None,
+                external_potential,
+                PoreSpecification::ChemicalPotential,
+            )?;
             profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));
 
             old_density = if let Some(Ok(l)) = profiles.last() {
@@ -245,7 +258,7 @@ where
 
     /// Calculate the phase transition from an empty to a filled pore.
     #[expect(clippy::too_many_arguments)]
-    pub fn phase_equilibrium<S: PoreSpecification<D>, X: Composition<f64, Dyn> + Clone>(
+    pub fn phase_equilibrium<S: Pore<D>, X: Composition<f64, Dyn> + Clone>(
         functional: &F,
         temperature: Temperature,
         p_min: Pressure,
@@ -261,8 +274,17 @@ where
         let bulk_init = State::new_npt(functional, temperature, p_max, x.clone(), Some(Liquid))?;
         let liquid_bulk = State::new_npt(functional, temperature, p_max, x.clone(), Some(Vapor))?;
 
-        let mut vapor = pore.initialize(&vapor_bulk, None, None)?.solve(solver)?;
-        let mut liquid = pore.initialize(&bulk_init, None, None)?.solve(solver)?;
+        let mut vapor = pore
+            .initialize(
+                &vapor_bulk,
+                None,
+                None,
+                PoreSpecification::ChemicalPotential,
+            )?
+            .solve(solver)?;
+        let mut liquid = pore
+            .initialize(&bulk_init, None, None, PoreSpecification::ChemicalPotential)?
+            .solve(solver)?;
 
         // calculate initial value for bulk density
         let n_dp_drho_v = (vapor.profile.moles() * vapor_bulk.dp_drho(Contributions::Total)).sum();

diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs
@@ -1,20 +1,20 @@
-use crate::WeightFunctionInfo;
 use crate::adsorption::{ExternalPotential, FluidParameters};
 use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, HelmholtzEnergyFunctionalDyn, MoleculeShape};
 use crate::functional_contribution::FunctionalContribution;
 use crate::geometry::{Axis, Geometry, Grid};
 use crate::profile::{DFTProfile, MAX_POTENTIAL};
 use crate::solver::DFTSolver;
+use crate::{DFTSpecification, WeightFunctionInfo};
 use feos_core::{Contributions, FeosResult, ReferenceSystem, ResidualDyn, State, StateHD};
 use nalgebra::{DVector, dvector};
 use ndarray::prelude::*;
 use ndarray::{Axis as Axis_nd, RemoveAxis};
 use num_dual::linalg::LU;
 use num_dual::{Dual64, DualNum};
 use quantity::{
-    _Moles, _Pressure, Density, Dimensionless, Energy, KELVIN, Length, MolarEnergy, Quantity, RGAS,
-    Temperature, Volume,
+    _Moles, _Pressure, Density, Dimensionless, Energy, KELVIN, Length, MolarEnergy, Moles,
+    Quantity, RGAS, Temperature, Volume,
 };
 use rustdct::DctNum;
 use std::ops::Sub;
@@ -52,14 +52,25 @@ impl Pore1D {
     }
 }
 
+/// Different ways that the thermodynamic state of the fluid in the pore
+/// can be specified.
+#[derive(Clone)]
+pub enum PoreSpecification {
+    /// Specify the chemical potential (via the bulk state).
+    ChemicalPotential,
+    /// Specify the amount of moles of every component.
+    Moles(Moles<Array1<f64>>),
+}
+
 /// Trait for the generic implementation of adsorption applications.
-pub trait PoreSpecification<D: Dimension> {
+pub trait Pore<D: Dimension> {
     /// Initialize a new single pore.
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<F>,
         density: Option<&Density<Array<f64, D::Larger>>>,
         external_potential: Option<&Array<f64, D::Larger>>,
+        specification: PoreSpecification,
     ) -> FeosResult<PoreProfile<D, F>>;
 
     /// Return the pore volume using Helium at 298 K as reference.
@@ -68,7 +79,7 @@ pub trait PoreSpecification<D: Dimension> {
         D::Larger: Dimension<Smaller = D>,
     {
         let bulk = State::new_pure(&&Helium, 298.0 * KELVIN, Density::from_reduced(1.0))?;
-        let pore = self.initialize(&bulk, None, None)?;
+        let pore = self.initialize(&bulk, None, None, PoreSpecification::ChemicalPotential)?;
         let pot = Dimensionless::from_reduced(
             pore.profile
                 .external_potential
@@ -96,6 +107,27 @@ where
     D::Smaller: Dimension<Larger = D>,
     <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
+    pub fn new(
+        grid: Grid,
+        bulk: &State<F>,
+        external_potential: Option<Array<f64, D::Larger>>,
+        density: Option<&Density<Array<f64, D::Larger>>>,
+        specification: PoreSpecification,
+    ) -> Self {
+        let mut profile = DFTProfile::new(grid, bulk, external_potential, density, Some(1));
+
+        // fix the number of particles
+        if let PoreSpecification::Moles(moles) = specification {
+            profile.specification = DFTSpecification::Moles(moles.to_reduced())
+        }
+
+        Self {
+            profile,
+            grand_potential: None,
+            interfacial_tension: None,
+        }
+    }
+
     pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> {
         // Solve the profile
         self.profile.solve(solver, debug)?;
@@ -180,12 +212,13 @@ where
     }
 }
 
-impl PoreSpecification<Ix1> for Pore1D {
+impl Pore<Ix1> for Pore1D {
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<F>,
         density: Option<&Density<Array2<f64>>>,
         external_potential: Option<&Array2<f64>>,
+        specification: PoreSpecification,
     ) -> FeosResult<PoreProfile1D<F>> {
         let dft: &F = &bulk.eos;
         let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS);
@@ -223,11 +256,13 @@ impl PoreSpecification<Ix1> for Pore1D {
         // initialize grid
         let grid = Grid::new_1d(axis);
 
-        Ok(PoreProfile {
-            profile: DFTProfile::new(grid, bulk, Some(external_potential), density, Some(1)),
-            grand_potential: None,
-            interfacial_tension: None,
-        })
+        Ok(PoreProfile::new(
+            grid,
+            bulk,
+            Some(external_potential),
+            density,
+            specification,
+        ))
     }
 }
 

diff --git a/crates/feos-dft/src/adsorption/pore2d.rs b/crates/feos-dft/src/adsorption/pore2d.rs
@@ -1,5 +1,5 @@
-use super::{FluidParameters, PoreProfile, PoreSpecification};
-use crate::{Axis, DFTProfile, Grid, HelmholtzEnergyFunctional};
+use super::{FluidParameters, Pore, PoreProfile};
+use crate::{Axis, Grid, HelmholtzEnergyFunctional, adsorption::pore::PoreSpecification};
 use feos_core::{FeosResult, State};
 use ndarray::{Array3, Ix2};
 use quantity::{Angle, Density, Length};
@@ -22,22 +22,25 @@ impl Pore2D {
     }
 }
 
-impl PoreSpecification<Ix2> for Pore2D {
+impl Pore<Ix2> for Pore2D {
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<F>,
         density: Option<&Density<Array3<f64>>>,
         external_potential: Option<&Array3<f64>>,
+        specification: PoreSpecification,
     ) -> FeosResult<PoreProfile<Ix2, F>> {
         // generate grid
         let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None);
         let y = Axis::new_cartesian(self.n_grid[1], self.system_size[1], None);
         let grid = Grid::Periodical2(x, y, self.angle);
 
-        Ok(PoreProfile {
-            profile: DFTProfile::new(grid, bulk, external_potential.cloned(), density, Some(1)),
-            grand_potential: None,
-            interfacial_tension: None,
-        })
+        Ok(PoreProfile::new(
+            grid,
+            bulk,
+            external_potential.cloned(),
+            density,
+            specification,
+        ))
     }
 }
diff --git a/crates/feos-dft/src/adsorption/pore3d.rs b/crates/feos-dft/src/adsorption/pore3d.rs
@@ -1,8 +1,8 @@
-use super::pore::{PoreProfile, PoreSpecification};
-use crate::adsorption::FluidParameters;
+use super::pore::{PoreProfile, Pore};
+use crate::adsorption::{FluidParameters, PoreSpecification};
 use crate::functional::HelmholtzEnergyFunctional;
 use crate::geometry::{Axis, Grid};
-use crate::profile::{CUTOFF_RADIUS, DFTProfile, MAX_POTENTIAL};
+use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL};
 use feos_core::{FeosError, FeosResult, ReferenceSystem, State};
 use ndarray::Zip;
 use ndarray::prelude::*;
@@ -48,12 +48,13 @@ impl Pore3D {
 /// Density profile and properties of a 3D confined system.
 pub type PoreProfile3D<F> = PoreProfile<Ix3, F>;
 
-impl PoreSpecification<Ix3> for Pore3D {
+impl Pore<Ix3> for Pore3D {
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<F>,
         density: Option<&Density<Array4<f64>>>,
         external_potential: Option<&Array4<f64>>,
+        specification: PoreSpecification,
     ) -> FeosResult<PoreProfile3D<F>> {
         let dft: &F = &bulk.eos;
 
@@ -94,11 +95,13 @@ impl PoreSpecification<Ix3> for Pore3D {
         )?;
         let grid = Grid::Periodical3(x, y, z, self.angles.unwrap_or([90.0 * DEGREES; 3]));
 
-        Ok(PoreProfile {
-            profile: DFTProfile::new(grid, bulk, Some(external_potential), density, Some(1)),
-            grand_potential: None,
-            interfacial_tension: None,
-        })
+        Ok(PoreProfile::new(
+            grid,
+            bulk,
+            Some(external_potential),
+            density,
+            specification,
+        ))
     }
 }
 

diff --git a/crates/feos-dft/src/interface/mod.rs b/crates/feos-dft/src/interface/mod.rs
@@ -2,12 +2,11 @@
 use crate::functional::HelmholtzEnergyFunctional;
 use crate::geometry::{Axis, Grid};
 use crate::pdgt::PdgtFunctionalProperties;
-use crate::profile::{DFTProfile, DFTSpecifications};
+use crate::profile::DFTProfile;
 use crate::solver::DFTSolver;
 use feos_core::{Contributions, FeosError, FeosResult, PhaseEquilibrium, ReferenceSystem};
 use ndarray::{Array1, Array2, Axis as Axis_nd, Ix1, s};
 use quantity::{Area, Density, Length, Moles, SurfaceTension, Temperature};
-use std::sync::Arc;
 
 mod surface_tension_diagram;
 pub use surface_tension_diagram::SurfaceTensionDiagram;
@@ -95,9 +94,7 @@ impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
 
         // specify specification
         if fix_equimolar_surface {
-            profile.profile.specification = Arc::new(DFTSpecifications::total_moles_from_profile(
-                &profile.profile,
-            ));
+            profile.profile.fix_total_moles();
         }
 
         profile
@@ -145,9 +142,7 @@ impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
 
         // specify specification
         if fix_equimolar_surface {
-            profile.profile.specification = Arc::new(DFTSpecifications::total_moles_from_profile(
-                &profile.profile,
-            ));
+            profile.profile.fix_total_moles();
         }
 
         Ok(profile)

diff --git a/crates/feos-dft/src/lib.rs b/crates/feos-dft/src/lib.rs
@@ -19,6 +19,6 @@ pub use functional::{HelmholtzEnergyFunctional, HelmholtzEnergyFunctionalDyn, Mo
 pub use functional_contribution::FunctionalContribution;
 pub use geometry::{Axis, Geometry, Grid};
 pub use pdgt::PdgtFunctionalProperties;
-pub use profile::{DFTProfile, DFTSpecification, DFTSpecifications};
+pub use profile::{DFTProfile, DFTSpecification};
 pub use solver::{DFTSolver, DFTSolverLog};
 pub use weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape};