Add calculation of reference pore volumes

feos-org · prehner · Feb 16, 2022 · Feb 16, 2022 · Feb 16, 2022 · f2453e645246b521a763273f38957f0e26e02724
commit f2453e645246b521a763273f38957f0e26e02724
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,8 @@ 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]
+### Added
+- The pore volume using Helium at 298 K as reference is now directly accesible from `Pore1D` and `Pore3D`. [#12](https://github.com/feos-org/feos-core/pull/12)
 
 ## [0.1.1] - 2022-02-14
 ### Added

diff --git a/examples/FundamentalMeasureTheory.ipynb b/examples/FundamentalMeasureTheory.ipynb
diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs
@@ -55,7 +55,7 @@ where
         })
     }
 
-    fn equilibrium<D: Dimension, F: HelmholtzEnergyFunctional>(
+    fn equilibrium<D: Dimension, F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         equilibrium: &Adsorption<U, D, F>,
     ) -> EosResult<(QuantityArray1<U>, QuantityArray1<U>)>
@@ -111,12 +111,12 @@ pub type Adsorption1D<U, F> = Adsorption<U, Ix1, F>;
 /// Container structure for adsorption isotherms in 3D pores.
 pub type Adsorption3D<U, F> = Adsorption<U, Ix3, F>;
 
-impl<U: EosUnit, D: Dimension, F: HelmholtzEnergyFunctional> Adsorption<U, D, F>
+impl<U: EosUnit, D: Dimension, F: HelmholtzEnergyFunctional + FluidParameters> Adsorption<U, D, F>
 where
     QuantityScalar<U>: std::fmt::Display,
     D::Larger: Dimension<Smaller = D>,
 {
-    fn new<S: PoreSpecification<U, D, F>>(
+    fn new<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
         pore: &S,
         profiles: Vec<EosResult<PoreProfile<U, D, F>>>,
@@ -129,7 +129,7 @@ where
     }
 
     /// Calculate an adsorption isotherm (starting at low pressure)
-    pub fn adsorption_isotherm<S: PoreSpecification<U, D, F>>(
+    pub fn adsorption_isotherm<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &PressureSpecification<U>,
@@ -142,7 +142,7 @@ where
     }
 
     /// Calculate an desorption isotherm (starting at high pressure)
-    pub fn desorption_isotherm<S: PoreSpecification<U, D, F>>(
+    pub fn desorption_isotherm<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &PressureSpecification<U>,
@@ -162,7 +162,7 @@ where
     }
 
     /// Calculate an equilibrium isotherm
-    pub fn equilibrium_isotherm<S: PoreSpecification<U, D, F>>(
+    pub fn equilibrium_isotherm<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &PressureSpecification<U>,
@@ -241,7 +241,7 @@ where
         }
     }
 
-    fn isotherm<S: PoreSpecification<U, D, F>>(
+    fn isotherm<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
         temperature: QuantityScalar<U>,
         pressure: &QuantityArray1<U>,
@@ -291,7 +291,7 @@ where
     }
 
     /// Calculate the phase transition from an empty to a filled pore.
-    pub fn phase_equilibrium<S: PoreSpecification<U, D, F>>(
+    pub fn phase_equilibrium<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
         temperature: QuantityScalar<U>,
         p_min: QuantityScalar<U>,

diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs
@@ -1,15 +1,17 @@
 use crate::adsorption::{ExternalPotential, FluidParameters};
 use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, DFT};
+use crate::functional_contribution::FunctionalContribution;
 use crate::geometry::{Axis, AxisGeometry, Grid};
 use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
 use crate::solver::DFTSolver;
-use feos_core::{Contributions, EosResult, EosUnit, State};
+use feos_core::{Contributions, EosResult, EosUnit, State, StateBuilder};
 use ndarray::prelude::*;
 use ndarray::Axis as Axis_nd;
 use ndarray::Zip;
 use ndarray_stats::QuantileExt;
 use quantity::{QuantityArray2, QuantityScalar};
+use std::rc::Rc;
 
 const POTENTIAL_OFFSET: f64 = 2.0;
 const DEFAULT_GRID_POINTS: usize = 2048;
@@ -76,16 +78,36 @@ impl<U> Pore3D<U> {
 }
 
 /// Trait for the generic implementation of adsorption applications.
-pub trait PoreSpecification<U, D: Dimension, F> {
+pub trait PoreSpecification<U: EosUnit, D: Dimension> {
     /// Initialize a new single pore.
-    fn initialize(
+    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array<f64, D::Larger>>,
     ) -> EosResult<PoreProfile<U, D, F>>;
 
     /// Return the number of spatial dimensions of the pore.
     fn dimension(&self) -> i32;
+
+    /// Return the pore volume using Helium at 298 K as reference.
+    fn pore_volume(&self) -> EosResult<QuantityScalar<U>>
+    where
+        D::Larger: Dimension<Smaller = D>,
+    {
+        let bulk = StateBuilder::new(&Rc::new(Helium::new()))
+            .temperature(298.0 * U::reference_temperature())
+            .volume(U::reference_volume())
+            .build()?;
+        let pore = self.initialize(&bulk, None)?;
+        let pot = pore
+            .profile
+            .external_potential
+            .index_axis(Axis(0), 0)
+            .mapv(|v| (-v).exp())
+            * U::reference_temperature()
+            / U::reference_temperature();
+        Ok(pore.profile.integrate(&pot))
+    }
 }
 
 /// Density profile and properties of a confined system in arbitrary dimensions.
@@ -149,10 +171,8 @@ where
     }
 }
 
-impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecification<U, Ix1, F>
-    for Pore1D<U>
-{
-    fn initialize(
+impl<U: EosUnit> PoreSpecification<U, Ix1> for Pore1D<U> {
+    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array2<f64>>,
@@ -203,10 +223,8 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecificati
     }
 }
 
-impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecification<U, Ix3, F>
-    for Pore3D<U>
-{
-    fn initialize(
+impl<U: EosUnit> PoreSpecification<U, Ix3> for Pore3D<U> {
+    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array4<f64>>,
@@ -420,3 +438,47 @@ fn calculate_distance2(
         rx.powi(2) + ry.powi(2) + rz.powi(2)
     })
 }
+
+const EPSILON_HE: f64 = 10.9;
+const SIGMA_HE: f64 = 2.64;
+
+struct Helium {
+    epsilon: Array1<f64>,
+    sigma: Array1<f64>,
+}
+
+impl Helium {
+    fn new() -> DFT<Self> {
+        let epsilon = arr1(&[EPSILON_HE]);
+        let sigma = arr1(&[SIGMA_HE]);
+        DFT::new_homosegmented(Self { epsilon, sigma }, &Array1::ones(1))
+    }
+}
+
+impl HelmholtzEnergyFunctional for Helium {
+    fn contributions(&self) -> &[Box<dyn FunctionalContribution>] {
+        &[]
+    }
+
+    fn subset(&self, _: &[usize]) -> DFT<Self> {
+        Self::new()
+    }
+
+    fn compute_max_density(&self, _: &Array1<f64>) -> f64 {
+        1.0
+    }
+}
+
+impl FluidParameters for Helium {
+    fn epsilon_k_ff(&self) -> Array1<f64> {
+        self.epsilon.clone()
+    }
+
+    fn sigma_ff(&self) -> &Array1<f64> {
+        &self.sigma
+    }
+
+    fn m(&self) -> Array1<f64> {
+        arr1(&[1.0])
+    }
+}
diff --git a/src/python/adsorption/external_potential.rs b/src/python/adsorption/external_potential.rs
@@ -222,7 +222,7 @@ impl PyExternalPotential {
 /// Returns
 /// -------
 /// Geometry
-#[pyclass(name = "Geometry", unsendable)]
+#[pyclass(name = "Geometry")]
 #[derive(Clone)]
 pub struct PyGeometry(pub AxisGeometry);
 

diff --git a/src/python/adsorption/pore.rs b/src/python/adsorption/pore.rs
@@ -74,6 +74,12 @@ macro_rules! impl_pore {
                     external_potential.map(|e| e.to_owned_array()).as_ref(),
                 )?))
             }
+
+            /// The pore volume using Helium at 298 K as reference.
+            #[getter]
+            fn get_pore_volume(&self) -> PyResult<PySINumber> {
+                Ok(self.0.pore_volume()?.into())
+            }
         }
 
         #[pymethods]