From b56eec6fc3eaee95dc915bed99ada64618db4b1e Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 24 Jul 2026 10:19:40 +0200 Subject: [PATCH] Pass external potentials as quantities and always cap them consistently --- CHANGELOG.md | 1 + crates/feos-derive/src/dft.rs | 2 +- .../feos-dft/src/adsorption/fea_potential.rs | 4 +- crates/feos-dft/src/adsorption/mod.rs | 2 +- crates/feos-dft/src/adsorption/pore.rs | 24 ++++-------- crates/feos-dft/src/adsorption/pore2d.rs | 4 +- crates/feos-dft/src/adsorption/pore3d.rs | 32 +++------------- crates/feos-dft/src/profile/mod.rs | 38 ++++++++++++++----- .../src/solvation/pair_correlation.rs | 12 ++---- .../src/solvation/solvation_profile.rs | 24 ++---------- crates/feos/src/hard_sphere/dft.rs | 7 ++-- crates/feos/src/pcsaft/dft/mod.rs | 17 +++++---- crates/feos/src/pets/dft/mod.rs | 24 +++++++----- crates/feos/src/saftvrqmie/dft/mod.rs | 12 +++--- py-feos/src/dft/adsorption/pore.rs | 25 ++++++------ py-feos/src/dft/profile.rs | 6 +-- py-feos/src/dft/solvation.rs | 7 +--- 17 files changed, 106 insertions(+), 135 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 05c507355..fe08b8fe4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - Changed data type of initial temperatures or pressure for phase equilibrium calculations (`TemperatureOrPressure::Other`) from `D` to `f64`. [#369](https://github.com/feos-org/feos/pull/369) - Reworked DFT solution algorithms slightly for the cases in which additional specifications are given. [#371](https://github.com/feos-org/feos/pull/371) +- External potentials are passed and returned as quantities (energies) instead of reduced units. [#372](https://github.com/feos-org/feos/pull/372) ### Removed - Removed the `DFTSpecification` trait in favor of only using the `DFTSpecification` enum (renamed from `DFTSpecifications`). [#371](https://github.com/feos-org/feos/pull/371) diff --git a/crates/feos-derive/src/dft.rs b/crates/feos-derive/src/dft.rs index d7dcd8e63..e5dd7a343 100644 --- a/crates/feos-derive/src/dft.rs +++ b/crates/feos-derive/src/dft.rs @@ -139,7 +139,7 @@ fn impl_pair_potential( } Ok(quote! { impl feos_dft::solvation::PairPotential for #ident { - fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> ndarray::Array2 { + fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> quantity::Energy> { match self { #(#pair_potential,)* } diff --git a/crates/feos-dft/src/adsorption/fea_potential.rs b/crates/feos-dft/src/adsorption/fea_potential.rs index d05ed6a01..f67e583f1 100644 --- a/crates/feos-dft/src/adsorption/fea_potential.rs +++ b/crates/feos-dft/src/adsorption/fea_potential.rs @@ -1,6 +1,6 @@ use super::pore3d::{calculate_distance2, evaluate_lj_potential}; use crate::Geometry; -use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL}; +use crate::profile::CUTOFF_RADIUS; use feos_core::ReferenceSystem; use gauss_quad::GaussLegendre; use ndarray::{Array1, Array2, Zip}; @@ -131,7 +131,7 @@ pub fn calculate_fea_potential( ) / temperature }) .sum(); - potential_2d[[i1, i2]] = (-potential_sum.min(MAX_POTENTIAL)).exp(); + potential_2d[[i1, i2]] = (-potential_sum).exp(); } } *f = (potential_2d * &weights).sum(); diff --git a/crates/feos-dft/src/adsorption/mod.rs b/crates/feos-dft/src/adsorption/mod.rs index 463989f67..3d11c4275 100644 --- a/crates/feos-dft/src/adsorption/mod.rs +++ b/crates/feos-dft/src/adsorption/mod.rs @@ -214,7 +214,7 @@ where .initialize(&bulk, None, None, PoreSpecification::ChemicalPotential)? .solve(solver)? .profile; - let external_potential = Some(&profile.external_potential); + let external_potential = Some(&profile.external_potential()); let mut old_density = Some(&profile.density); for i in 0..pressure.len() { diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs index 5db29b1f7..493229289 100644 --- a/crates/feos-dft/src/adsorption/pore.rs +++ b/crates/feos-dft/src/adsorption/pore.rs @@ -3,7 +3,7 @@ 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::profile::DFTProfile; use crate::solver::DFTSolver; use crate::{DFTSpecification, WeightFunctionInfo}; use feos_core::{Contributions, FeosResult, ReferenceSystem, ResidualDyn, State, StateHD}; @@ -69,7 +69,7 @@ pub trait Pore { &self, bulk: &State, density: Option<&Density>>, - external_potential: Option<&Array>, + external_potential: Option<&Energy>>, specification: PoreSpecification, ) -> FeosResult>; @@ -110,7 +110,7 @@ where pub fn new( grid: Grid, bulk: &State, - external_potential: Option>, + external_potential: Option>>, density: Option<&Density>>, specification: PoreSpecification, ) -> Self { @@ -217,7 +217,7 @@ impl Pore for Pore1D { &self, bulk: &State, density: Option<&Density>>, - external_potential: Option<&Array2>, + external_potential: Option<&Energy>>, specification: PoreSpecification, ) -> FeosResult> { let dft: &F = &bulk.eos; @@ -247,7 +247,6 @@ impl Pore for Pore1D { &self.potential, dft, &axis, - self.potential_cutoff, ) }, |e| e.clone(), @@ -272,9 +271,7 @@ fn external_potential_1d( potential: &ExternalPotential, fluid_parameters: &P, axis: &Axis, - potential_cutoff: Option, -) -> Array2 { - let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL); +) -> Energy> { let effective_pore_size = match axis.geometry { Geometry::Spherical => pore_width.to_reduced(), Geometry::Cylindrical => pore_width.to_reduced(), @@ -305,21 +302,16 @@ fn external_potential_1d( fluid_parameters, t, ), - } / t; + }; for (i, &z) in axis.grid.iter().enumerate() { if z > effective_pore_size { external_potential .index_axis_mut(Axis_nd(1), i) - .fill(potential_cutoff); + .fill(f64::INFINITY); } } - external_potential.map_inplace(|x| { - if *x > potential_cutoff { - *x = potential_cutoff - } - }); - external_potential + Energy::from_reduced(external_potential) } const EPSILON_HE: f64 = 10.9; diff --git a/crates/feos-dft/src/adsorption/pore2d.rs b/crates/feos-dft/src/adsorption/pore2d.rs index f7290c6e3..396434187 100644 --- a/crates/feos-dft/src/adsorption/pore2d.rs +++ b/crates/feos-dft/src/adsorption/pore2d.rs @@ -2,7 +2,7 @@ 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}; +use quantity::{Angle, Density, Energy, Length}; pub struct Pore2D { system_size: [Length; 2], @@ -27,7 +27,7 @@ impl Pore for Pore2D { &self, bulk: &State, density: Option<&Density>>, - external_potential: Option<&Array3>, + external_potential: Option<&Energy>>, specification: PoreSpecification, ) -> FeosResult> { // generate grid diff --git a/crates/feos-dft/src/adsorption/pore3d.rs b/crates/feos-dft/src/adsorption/pore3d.rs index bfbfd38df..fdf02800a 100644 --- a/crates/feos-dft/src/adsorption/pore3d.rs +++ b/crates/feos-dft/src/adsorption/pore3d.rs @@ -1,12 +1,12 @@ -use super::pore::{PoreProfile, Pore}; +use super::pore::{Pore, PoreProfile}; use crate::adsorption::{FluidParameters, PoreSpecification}; use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::{Axis, Grid}; -use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL}; +use crate::profile::CUTOFF_RADIUS; use feos_core::{FeosError, FeosResult, ReferenceSystem, State}; use ndarray::Zip; use ndarray::prelude::*; -use quantity::{Angle, DEGREES, Density, Length}; +use quantity::{Angle, DEGREES, Density, Energy, Length}; /// Parameters required to specify a 3D pore. pub struct Pore3D { @@ -16,12 +16,10 @@ pub struct Pore3D { coordinates: Length>, sigma_ss: Array1, epsilon_k_ss: Array1, - potential_cutoff: Option, cutoff_radius: Option, } impl Pore3D { - #[expect(clippy::too_many_arguments)] pub fn new( system_size: [Length; 3], n_grid: [usize; 3], @@ -29,7 +27,6 @@ impl Pore3D { sigma_ss: Array1, epsilon_k_ss: Array1, angles: Option<[Angle; 3]>, - potential_cutoff: Option, cutoff_radius: Option, ) -> Self { Self { @@ -39,7 +36,6 @@ impl Pore3D { coordinates, sigma_ss, epsilon_k_ss, - potential_cutoff, cutoff_radius, } } @@ -53,7 +49,7 @@ impl Pore for Pore3D { &self, bulk: &State, density: Option<&Density>>, - external_potential: Option<&Array4>, + external_potential: Option<&Energy>>, specification: PoreSpecification, ) -> FeosResult> { let dft: &F = &bulk.eos; @@ -65,9 +61,6 @@ impl Pore for Pore3D { let coordinates = self.coordinates.to_reduced(); - // temperature - let t = bulk.temperature.to_reduced(); - // For non-orthorombic unit cells, the external potential has to be // provided at the moment if let (Some(_), None) = (self.angles, external_potential) { @@ -87,8 +80,6 @@ impl Pore for Pore3D { &self.sigma_ss, &self.epsilon_k_ss, self.cutoff_radius, - self.potential_cutoff, - t, ) }, |e| Ok(e.clone()), @@ -105,7 +96,6 @@ impl Pore for Pore3D { } } -#[expect(clippy::too_many_arguments)] pub fn external_potential_3d( functional: &F, axis: [&Axis; 3], @@ -114,9 +104,7 @@ pub fn external_potential_3d( sigma_ss: &Array1, epsilon_ss: &Array1, cutoff_radius: Option, - potential_cutoff: Option, - reduced_temperature: f64, -) -> FeosResult> { +) -> FeosResult>> { // allocate external potential let m = functional.m(); let mut external_potential = Array4::zeros(( @@ -167,17 +155,9 @@ pub fn external_potential_3d( ) }) .sum::() - / reduced_temperature - }); - - let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL); - external_potential.map_inplace(|x| { - if *x > potential_cutoff { - *x = potential_cutoff - } }); - Ok(external_potential) + Ok(Energy::from_reduced(external_potential)) } /// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index 26b1d0d94..c9423fded 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -9,13 +9,15 @@ use ndarray::{ RemoveAxis, }; use num_dual::DualNum; -use quantity::{_Volume, DEGREES, Density, Length, Moles, Quantity, Temperature, Volume}; +use quantity::{ + _Volume, DEGREES, Density, Energy, Entropy, Length, Moles, Quantity, Temperature, Volume, +}; use std::ops::{Add, MulAssign}; use std::sync::Arc; mod properties; -pub(crate) const MAX_POTENTIAL: f64 = 50.0; +const MAX_POTENTIAL: f64 = 50.0; #[cfg(feature = "rayon")] pub(crate) const CUTOFF_RADIUS: f64 = 14.0; @@ -161,7 +163,7 @@ where pub fn new( grid: Grid, bulk: &State, - external_potential: Option>, + external_potential: Option>>, density: Option<&Density>>, lanczos: Option, ) -> Self { @@ -171,13 +173,24 @@ where let convolver = ConvolverFFT::plan(&grid, &weight_functions, lanczos); // initialize external potential - let external_potential = external_potential.unwrap_or_else(|| { - let mut n_grid = vec![bulk.eos.component_index().len()]; - grid.axes() - .iter() - .for_each(|&ax| n_grid.push(ax.grid.len())); - Array::zeros(n_grid).into_dimensionality().unwrap() - }); + let external_potential = external_potential.map_or_else( + || { + let mut n_grid = vec![bulk.eos.component_index().len()]; + grid.axes() + .iter() + .for_each(|&ax| n_grid.push(ax.grid.len())); + Array::zeros(n_grid).into_dimensionality().unwrap() + }, + |e| { + let mut external_potential = e.into_reduced() / t; + external_potential.map_inplace(|x| { + if *x > MAX_POTENTIAL { + *x = MAX_POTENTIAL + } + }); + external_potential + }, + ); // initialize density let density = if let Some(density) = density { @@ -223,6 +236,11 @@ where let moles = self.integrate_reduced_comp(&rho).sum(); self.specification = DFTSpecification::TotalMoles(moles); } + + /// Return the external potential in SI units. + pub fn external_potential(&self) -> Energy> { + Entropy::from_reduced(self.external_potential.clone()) * self.temperature + } } impl DFTProfile diff --git a/crates/feos-dft/src/solvation/pair_correlation.rs b/crates/feos-dft/src/solvation/pair_correlation.rs index 25849af4f..85b435d2e 100644 --- a/crates/feos-dft/src/solvation/pair_correlation.rs +++ b/crates/feos-dft/src/solvation/pair_correlation.rs @@ -1,6 +1,5 @@ //! Functionalities for the calculation of pair correlation functions. use crate::functional::HelmholtzEnergyFunctional; -use crate::profile::MAX_POTENTIAL; use crate::solver::DFTSolver; use crate::{Axis, DFTProfile, Grid}; use feos_core::{Contributions, FeosResult, ReferenceSystem, State}; @@ -12,11 +11,11 @@ use std::ops::Deref; /// models. pub trait PairPotential { /// Return the pair potential of particle i with all other particles. - fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Array2; + fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Energy>; } impl, T: PairPotential> PairPotential for C { - fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Array2 { + fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Energy> { T::pair_potential(self, i, r, temperature) } } @@ -38,12 +37,7 @@ impl PairCorrelation { // calculate external potential let t = bulk.temperature.to_reduced(); - let mut external_potential = dft.pair_potential(test_particle, &axis.grid, t) / t; - external_potential.map_inplace(|x| { - if *x > MAX_POTENTIAL { - *x = MAX_POTENTIAL - } - }); + let external_potential = dft.pair_potential(test_particle, &axis.grid, t) / t; let grid = Grid::Spherical(axis); Self { diff --git a/crates/feos-dft/src/solvation/solvation_profile.rs b/crates/feos-dft/src/solvation/solvation_profile.rs index 8785876d5..e31ca6494 100644 --- a/crates/feos-dft/src/solvation/solvation_profile.rs +++ b/crates/feos-dft/src/solvation/solvation_profile.rs @@ -1,7 +1,7 @@ use crate::adsorption::FluidParameters; use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::{Axis, Grid}; -use crate::profile::{CUTOFF_RADIUS, DFTProfile, MAX_POTENTIAL}; +use crate::profile::{CUTOFF_RADIUS, DFTProfile}; use crate::solver::DFTSolver; use feos_core::{Contributions, FeosResult, ReferenceSystem, State}; use ndarray::Zip; @@ -40,7 +40,6 @@ impl SolvationProfile { } impl SolvationProfile { - #[expect(clippy::too_many_arguments)] pub fn new( bulk: &State, n_grid: [usize; 3], @@ -49,7 +48,6 @@ impl SolvationProfile { epsilon_ss: Array1, system_size: Option<[Length; 3]>, cutoff_radius: Option, - potential_cutoff: Option, ) -> FeosResult { let dft: &F = &bulk.eos; @@ -77,9 +75,6 @@ impl SolvationProfile { coordinates = coordinates + shift; - // temperature - let t = bulk.temperature.to_reduced(); - // calculate external potential let external_potential = external_potential_3d( dft, @@ -88,8 +83,6 @@ impl SolvationProfile { sigma_ss, epsilon_ss, cutoff_radius, - potential_cutoff, - t, )?; let grid = Grid::Cartesian3(x, y, z); @@ -102,7 +95,6 @@ impl SolvationProfile { } } -#[expect(clippy::too_many_arguments)] fn external_potential_3d( functional: &F, axis: [&Axis; 3], @@ -110,9 +102,7 @@ fn external_potential_3d( sigma_ss: Array1, epsilon_ss: Array1, cutoff_radius: Option, - potential_cutoff: Option, - reduced_temperature: f64, -) -> FeosResult> { +) -> FeosResult>> { // allocate external potential let m = functional.m(); let mut external_potential = Array4::zeros(( @@ -150,17 +140,9 @@ fn external_potential_3d( ) }) .sum::() - / reduced_temperature - }); - - let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL); - external_potential.map_inplace(|x| { - if *x > potential_cutoff { - *x = potential_cutoff - } }); - Ok(external_potential) + Ok(Energy::from_reduced(external_potential)) } /// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment diff --git a/crates/feos/src/hard_sphere/dft.rs b/crates/feos/src/hard_sphere/dft.rs index ba8edf25a..d48d6327a 100644 --- a/crates/feos/src/hard_sphere/dft.rs +++ b/crates/feos/src/hard_sphere/dft.rs @@ -8,6 +8,7 @@ use feos_dft::{ use nalgebra::DVector; use ndarray::*; use num_dual::DualNum; +use quantity::Energy; use std::f64::consts::PI; use super::{HardSphereProperties, MonomerShape}; @@ -352,15 +353,15 @@ impl HelmholtzEnergyFunctionalDyn for FMTFunctional { } impl PairPotential for FMTFunctional { - fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Array2 { + fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Energy> { let s = &self.properties.sigma; - Array::from_shape_fn((s.len(), r.len()), |(j, k)| { + Energy::new(Array::from_shape_fn((s.len(), r.len()), |(j, k)| { if r[k] > 0.5 * (s[i] + s[j]) { 0.0 } else { f64::INFINITY } - }) + })) } } diff --git a/crates/feos/src/pcsaft/dft/mod.rs b/crates/feos/src/pcsaft/dft/mod.rs index b3afa5435..52c16aa4a 100644 --- a/crates/feos/src/pcsaft/dft/mod.rs +++ b/crates/feos/src/pcsaft/dft/mod.rs @@ -3,7 +3,7 @@ use super::parameters::PcSaftPars; use crate::association::{Association, YuWuAssociationFunctional}; use crate::hard_sphere::{FMTContribution, FMTVersion}; use crate::pcsaft::eos::PcSaftOptions; -use feos_core::{FeosResult, Molarweight, ResidualDyn, StateHD, Subset}; +use feos_core::{FeosResult, Molarweight, ReferenceSystem, ResidualDyn, StateHD, Subset}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; @@ -14,7 +14,7 @@ use nalgebra::DVector; use ndarray::{Array1, Array2}; use num_dual::DualNum; use num_traits::One; -use quantity::MolarWeight; +use quantity::{Energy, MolarWeight}; use std::f64::consts::FRAC_PI_6; mod dispersion; @@ -164,13 +164,16 @@ impl FluidParameters for PcSaftFunctional { } impl PairPotential for PcSaftFunctional { - fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Array2 { + fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Energy> { let sigma_ij = &self.params.sigma_ij; let eps_ij_4 = 4.0 * &self.params.epsilon_k_ij; - Array2::from_shape_fn((self.params.m.len(), r.len()), |(j, k)| { - let att = (sigma_ij[(i, j)] / r[k]).powi(6); - eps_ij_4[(i, j)] * att * (att - 1.0) - }) + Energy::from_reduced(Array2::from_shape_fn( + (self.params.m.len(), r.len()), + |(j, k)| { + let att = (sigma_ij[(i, j)] / r[k]).powi(6); + eps_ij_4[(i, j)] * att * (att - 1.0) + }, + )) } } diff --git a/crates/feos/src/pets/dft/mod.rs b/crates/feos/src/pets/dft/mod.rs index 14215fee9..ef7e5d799 100644 --- a/crates/feos/src/pets/dft/mod.rs +++ b/crates/feos/src/pets/dft/mod.rs @@ -3,7 +3,7 @@ use super::eos::PetsOptions; use super::parameters::PetsParameters; use crate::hard_sphere::{FMTContribution, FMTVersion}; use dispersion::AttractiveFunctional; -use feos_core::FeosResult; +use feos_core::{FeosResult, ReferenceSystem}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; @@ -12,6 +12,7 @@ use nalgebra::DVector; use ndarray::{Array1, Array2}; use num_dual::DualNum; use pure_pets_functional::*; +use quantity::Energy; mod dispersion; mod pure_pets_functional; @@ -75,18 +76,21 @@ impl FluidParameters for Pets { } impl PairPotential for Pets { - fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Array2 { + fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Energy> { let eps_ij_4 = 4.0 * self.epsilon_k_ij.clone(); let shift_ij = &eps_ij_4 * (2.5.powi(-12) - 2.5.powi(-6)); let rc_ij = 2.5 * &self.sigma_ij; - Array2::from_shape_fn((self.sigma.len(), r.len()), |(j, k)| { - if r[k] > rc_ij[(i, j)] { - 0.0 - } else { - let att = (self.sigma_ij[(i, j)] / r[k]).powi(6); - eps_ij_4[(i, j)] * att * (att - 1.0) - shift_ij[(i, j)] - } - }) + Energy::from_reduced(Array2::from_shape_fn( + (self.sigma.len(), r.len()), + |(j, k)| { + if r[k] > rc_ij[(i, j)] { + 0.0 + } else { + let att = (self.sigma_ij[(i, j)] / r[k]).powi(6); + eps_ij_4[(i, j)] * att * (att - 1.0) - shift_ij[(i, j)] + } + }, + )) } } diff --git a/crates/feos/src/saftvrqmie/dft/mod.rs b/crates/feos/src/saftvrqmie/dft/mod.rs index 8f29059b5..6cde777b4 100644 --- a/crates/feos/src/saftvrqmie/dft/mod.rs +++ b/crates/feos/src/saftvrqmie/dft/mod.rs @@ -3,7 +3,7 @@ use crate::hard_sphere::{FMTContribution, FMTVersion, HardSphereProperties, Mono use crate::saftvrqmie::eos::SaftVRQMieOptions; use crate::saftvrqmie::parameters::{SaftVRQMieParameters, SaftVRQMiePars}; use dispersion::AttractiveFunctional; -use feos_core::FeosResult; +use feos_core::{FeosResult, ReferenceSystem}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; @@ -12,6 +12,7 @@ use nalgebra::DVector; use ndarray::{Array, Array1, Array2}; use non_additive_hs::NonAddHardSphereFunctional; use num_dual::DualNum; +use quantity::Energy; mod dispersion; mod non_additive_hs; @@ -79,10 +80,11 @@ impl FluidParameters for SaftVRQMie { } impl PairPotential for SaftVRQMie { - fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Array2 { - Array::from_shape_fn((self.params.m.len(), r.len()), |(j, k)| { - self.params.qmie_potential_ij(i, j, r[k], temperature)[0] - }) + fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Energy> { + Energy::from_reduced(Array::from_shape_fn( + (self.params.m.len(), r.len()), + |(j, k)| self.params.qmie_potential_ij(i, j, r[k], temperature)[0], + )) } } diff --git a/py-feos/src/dft/adsorption/pore.rs b/py-feos/src/dft/adsorption/pore.rs index 7668bebfb..6e39d57ab 100644 --- a/py-feos/src/dft/adsorption/pore.rs +++ b/py-feos/src/dft/adsorption/pore.rs @@ -139,7 +139,7 @@ impl PyPore1D { /// The bulk state in equilibrium with the pore. /// density : SIArray2, optional /// Initial values for the density profile. - /// external_potential : numpy.ndarray[float], optional + /// external_potential : SIArray2, optional /// The external potential in the pore. Used to /// save computation time in the case of costly /// evaluations of external potentials. @@ -158,7 +158,7 @@ impl PyPore1D { &self, bulk: &PyState, density: Option>>, - external_potential: Option<&Bound<'_, PyArray2>>, + external_potential: Option>>, specification: PyPoreSpecification, ) -> PyResult { Ok(PyPoreProfile1D( @@ -166,7 +166,7 @@ impl PyPore1D { .initialize( &bulk.0, density.as_ref(), - external_potential.map(|e| e.to_owned_array()).as_ref(), + external_potential.as_ref(), specification.0, ) .map_err(PyFeosError::from)?, @@ -236,7 +236,7 @@ impl PyPore2D { /// The bulk state in equilibrium with the pore. /// density : SIArray3, optional /// Initial values for the density profile. - /// external_potential : numpy.ndarray[float], optional + /// external_potential : SIArray3, optional /// The external potential in the pore. Used to /// save computation time in the case of costly /// evaluations of external potentials. @@ -255,7 +255,7 @@ impl PyPore2D { &self, bulk: &PyState, density: Option>>, - external_potential: Option<&Bound<'_, PyArray3>>, + external_potential: Option>>, specification: PyPoreSpecification, ) -> PyResult { Ok(PyPoreProfile2D( @@ -263,7 +263,7 @@ impl PyPore2D { .initialize( &bulk.0, density.as_ref(), - external_potential.map(|e| e.to_owned_array()).as_ref(), + external_potential.as_ref(), specification.0, ) .map_err(PyFeosError::from)?, @@ -325,10 +325,9 @@ impl_pore_profile!(PyPoreProfile3D); impl PyPore3D { #[new] #[pyo3( - text_signature = "(system_size, n_grid, coordinates, sigma_ss, epsilon_k_ss, angles=None, potential_cutoff=None, cutoff_radius=None)" + text_signature = "(system_size, n_grid, coordinates, sigma_ss, epsilon_k_ss, angles=None, cutoff_radius=None)" )] - #[pyo3(signature = (system_size, n_grid, coordinates, sigma_ss, epsilon_k_ss, angles=None, potential_cutoff=None, cutoff_radius=None))] - #[expect(clippy::too_many_arguments)] + #[pyo3(signature = (system_size, n_grid, coordinates, sigma_ss, epsilon_k_ss, angles=None, cutoff_radius=None))] fn new( system_size: [Length; 3], n_grid: [usize; 3], @@ -336,7 +335,6 @@ impl PyPore3D { sigma_ss: &Bound<'_, PyArray1>, epsilon_k_ss: &Bound<'_, PyArray1>, angles: Option<[Angle; 3]>, - potential_cutoff: Option, cutoff_radius: Option, ) -> Self { Self(Pore3D::new( @@ -346,7 +344,6 @@ impl PyPore3D { sigma_ss.to_owned_array(), epsilon_k_ss.to_owned_array(), angles, - potential_cutoff, cutoff_radius, )) } @@ -359,7 +356,7 @@ impl PyPore3D { /// The bulk state in equilibrium with the pore. /// density : SIArray4, optional /// Initial values for the density profile. - /// external_potential : numpy.ndarray[float], optional + /// external_potential : SIArray4, optional /// The external potential in the pore. Used to /// save computation time in the case of costly /// evaluations of external potentials. @@ -378,7 +375,7 @@ impl PyPore3D { &self, bulk: &PyState, density: Option>>, - external_potential: Option<&Bound<'_, PyArray4>>, + external_potential: Option>>, specification: PyPoreSpecification, ) -> PyResult { Ok(PyPoreProfile3D( @@ -386,7 +383,7 @@ impl PyPore3D { .initialize( &bulk.0, density.as_ref(), - external_potential.map(|e| e.to_owned_array()).as_ref(), + external_potential.as_ref(), specification.0, ) .map_err(PyFeosError::from)?, diff --git a/py-feos/src/dft/profile.rs b/py-feos/src/dft/profile.rs index 68c4aec9c..b60fcafdb 100644 --- a/py-feos/src/dft/profile.rs +++ b/py-feos/src/dft/profile.rs @@ -74,8 +74,8 @@ macro_rules! impl_profile { } #[getter] - fn get_external_potential<'py>(&self, py: Python<'py>) -> Bound<'py, $py_arr2> { - self.0.profile.external_potential.view().to_pyarray(py) + fn get_external_potential(&self) -> Energy<$si_arr2> { + self.0.profile.external_potential().clone() } #[getter] @@ -272,6 +272,6 @@ macro_rules! impl_3d_profile { }; } -pub(crate) use {impl_1d_profile, impl_2d_profile, impl_profile}; #[cfg(feature = "rayon")] pub(crate) use impl_3d_profile; +pub(crate) use {impl_1d_profile, impl_2d_profile, impl_profile}; diff --git a/py-feos/src/dft/solvation.rs b/py-feos/src/dft/solvation.rs index 4dc012972..9509d6b99 100644 --- a/py-feos/src/dft/solvation.rs +++ b/py-feos/src/dft/solvation.rs @@ -55,10 +55,9 @@ impl_3d_profile!(PySolvationProfile, get_x, get_y, get_z); impl PySolvationProfile { #[new] #[pyo3( - text_signature = "(bulk, n_grid, coordinates, sigma, epsilon_k, system_size=None, cutoff_radius=None, potential_cutoff=None)" + text_signature = "(bulk, n_grid, coordinates, sigma, epsilon_k, system_size=None, cutoff_radius=None)" )] - #[pyo3(signature = (bulk, n_grid, coordinates, sigma, epsilon_k, system_size=None, cutoff_radius=None, potential_cutoff=None))] - #[expect(clippy::too_many_arguments)] + #[pyo3(signature = (bulk, n_grid, coordinates, sigma, epsilon_k, system_size=None, cutoff_radius=None))] fn new<'py>( bulk: &PyState, n_grid: [usize; 3], @@ -67,7 +66,6 @@ impl PySolvationProfile { epsilon_k: &Bound<'py, PyArray1>, system_size: Option<[Length; 3]>, cutoff_radius: Option, - potential_cutoff: Option, ) -> PyResult { Ok(Self( SolvationProfile::new( @@ -78,7 +76,6 @@ impl PySolvationProfile { epsilon_k.to_owned_array(), system_size, cutoff_radius, - potential_cutoff, ) .map_err(PyFeosError::from)?, ))