Skip to content

Pass external potentials as quantities and always cap them consistently #372

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prehner merged 1 commit into from
Jul 24, 2026
Merged

Pass external potentials as quantities and always cap them consistently #372

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion crates/feos-derive/src/dft.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<f64>, temperature: f64) -> ndarray::Array2<f64> {
fn pair_potential(&self, i: usize, r: &Array1<f64>, temperature: f64) -> quantity::Energy<ndarray::Array2<f64>> {
match self {
#(#pair_potential,)*
}
Expand Down
4 changes: 2 additions & 2 deletions crates/feos-dft/src/adsorption/fea_potential.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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};
Expand Down Expand Up @@ -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();
Expand Down
2 changes: 1 addition & 1 deletion crates/feos-dft/src/adsorption/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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() {
Expand Down
24 changes: 8 additions & 16 deletions crates/feos-dft/src/adsorption/pore.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -69,7 +69,7 @@ pub trait Pore<D: Dimension> {
&self,
bulk: &State<F>,
density: Option<&Density<Array<f64, D::Larger>>>,
external_potential: Option<&Array<f64, D::Larger>>,
external_potential: Option<&Energy<Array<f64, D::Larger>>>,
specification: PoreSpecification,
) -> FeosResult<PoreProfile<D, F>>;

Expand Down Expand Up @@ -110,7 +110,7 @@ where
pub fn new(
grid: Grid,
bulk: &State<F>,
external_potential: Option<Array<f64, D::Larger>>,
external_potential: Option<Energy<Array<f64, D::Larger>>>,
density: Option<&Density<Array<f64, D::Larger>>>,
specification: PoreSpecification,
) -> Self {
Expand Down Expand Up @@ -217,7 +217,7 @@ impl Pore<Ix1> for Pore1D {
&self,
bulk: &State<F>,
density: Option<&Density<Array2<f64>>>,
external_potential: Option<&Array2<f64>>,
external_potential: Option<&Energy<Array2<f64>>>,
specification: PoreSpecification,
) -> FeosResult<PoreProfile1D<F>> {
let dft: &F = &bulk.eos;
Expand Down Expand Up @@ -247,7 +247,6 @@ impl Pore<Ix1> for Pore1D {
&self.potential,
dft,
&axis,
self.potential_cutoff,
)
},
|e| e.clone(),
Expand All @@ -272,9 +271,7 @@ fn external_potential_1d<P: HelmholtzEnergyFunctional + FluidParameters>(
potential: &ExternalPotential,
fluid_parameters: &P,
axis: &Axis,
potential_cutoff: Option<f64>,
) -> Array2<f64> {
let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL);
) -> Energy<Array2<f64>> {
let effective_pore_size = match axis.geometry {
Geometry::Spherical => pore_width.to_reduced(),
Geometry::Cylindrical => pore_width.to_reduced(),
Expand Down Expand Up @@ -305,21 +302,16 @@ fn external_potential_1d<P: HelmholtzEnergyFunctional + FluidParameters>(
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;
Expand Down
4 changes: 2 additions & 2 deletions crates/feos-dft/src/adsorption/pore2d.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<f64>; 2],
Expand All @@ -27,7 +27,7 @@ impl Pore<Ix2> for Pore2D {
&self,
bulk: &State<F>,
density: Option<&Density<Array3<f64>>>,
external_potential: Option<&Array3<f64>>,
external_potential: Option<&Energy<Array3<f64>>>,
specification: PoreSpecification,
) -> FeosResult<PoreProfile<Ix2, F>> {
// generate grid
Expand Down
32 changes: 6 additions & 26 deletions crates/feos-dft/src/adsorption/pore3d.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 {
Expand All @@ -16,20 +16,17 @@ pub struct Pore3D {
coordinates: Length<Array2<f64>>,
sigma_ss: Array1<f64>,
epsilon_k_ss: Array1<f64>,
potential_cutoff: Option<f64>,
cutoff_radius: Option<Length>,
}

impl Pore3D {
#[expect(clippy::too_many_arguments)]
pub fn new(
system_size: [Length; 3],
n_grid: [usize; 3],
coordinates: Length<Array2<f64>>,
sigma_ss: Array1<f64>,
epsilon_k_ss: Array1<f64>,
angles: Option<[Angle; 3]>,
potential_cutoff: Option<f64>,
cutoff_radius: Option<Length>,
) -> Self {
Self {
Expand All @@ -39,7 +36,6 @@ impl Pore3D {
coordinates,
sigma_ss,
epsilon_k_ss,
potential_cutoff,
cutoff_radius,
}
}
Expand All @@ -53,7 +49,7 @@ impl Pore<Ix3> for Pore3D {
&self,
bulk: &State<F>,
density: Option<&Density<Array4<f64>>>,
external_potential: Option<&Array4<f64>>,
external_potential: Option<&Energy<Array4<f64>>>,
specification: PoreSpecification,
) -> FeosResult<PoreProfile3D<F>> {
let dft: &F = &bulk.eos;
Expand All @@ -65,9 +61,6 @@ impl Pore<Ix3> 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) {
Expand All @@ -87,8 +80,6 @@ impl Pore<Ix3> for Pore3D {
&self.sigma_ss,
&self.epsilon_k_ss,
self.cutoff_radius,
self.potential_cutoff,
t,
)
},
|e| Ok(e.clone()),
Expand All @@ -105,7 +96,6 @@ impl Pore<Ix3> for Pore3D {
}
}

#[expect(clippy::too_many_arguments)]
pub fn external_potential_3d<F: HelmholtzEnergyFunctional + FluidParameters>(
functional: &F,
axis: [&Axis; 3],
Expand All @@ -114,9 +104,7 @@ pub fn external_potential_3d<F: HelmholtzEnergyFunctional + FluidParameters>(
sigma_ss: &Array1<f64>,
epsilon_ss: &Array1<f64>,
cutoff_radius: Option<Length>,
potential_cutoff: Option<f64>,
reduced_temperature: f64,
) -> FeosResult<Array4<f64>> {
) -> FeosResult<Energy<Array4<f64>>> {
// allocate external potential
let m = functional.m();
let mut external_potential = Array4::zeros((
Expand Down Expand Up @@ -167,17 +155,9 @@ pub fn external_potential_3d<F: HelmholtzEnergyFunctional + FluidParameters>(
)
})
.sum::<f64>()
/ 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
Expand Down
38 changes: 28 additions & 10 deletions crates/feos-dft/src/profile/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -161,7 +163,7 @@ where
pub fn new(
grid: Grid,
bulk: &State<F>,
external_potential: Option<Array<f64, D::Larger>>,
external_potential: Option<Energy<Array<f64, D::Larger>>>,
density: Option<&Density<Array<f64, D::Larger>>>,
lanczos: Option<i32>,
) -> Self {
Expand All @@ -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 {
Expand Down Expand Up @@ -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<Array<f64, D::Larger>> {
Entropy::from_reduced(self.external_potential.clone()) * self.temperature
}
}

impl<D: Dimension, F: HelmholtzEnergyFunctional> DFTProfile<D, F>
Expand Down
12 changes: 3 additions & 9 deletions crates/feos-dft/src/solvation/pair_correlation.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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};
Expand All @@ -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<f64>, temperature: f64) -> Array2<f64>;
fn pair_potential(&self, i: usize, r: &Array1<f64>, temperature: f64) -> Energy<Array2<f64>>;
}

impl<C: Deref<Target = T>, T: PairPotential> PairPotential for C {
fn pair_potential(&self, i: usize, r: &Array1<f64>, temperature: f64) -> Array2<f64> {
fn pair_potential(&self, i: usize, r: &Array1<f64>, temperature: f64) -> Energy<Array2<f64>> {
T::pair_potential(self, i, r, temperature)
}
}
Expand All @@ -38,12 +37,7 @@ impl<F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<F> {

// 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 {
Expand Down
Loading