Ideal gas properties from DFT

feos-org · prehner · Dec 10, 2024 · Dec 5, 2024 · Jan 25, 2024 · Jan 25, 2024
commit a6dfc8d7018e43a28c591e93970fcb18565ce33a
diff --git a/Cargo.toml b/Cargo.toml
@@ -70,7 +70,7 @@ codegen-units = 1
 
 
 [features]
-default = []
+default = ["python", "pcsaft", "dft"]
 dft = ["feos-dft", "petgraph"]
 estimator = []
 association = []

diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs
@@ -6,19 +6,27 @@ use crate::geometry::{Axis, Geometry, Grid};
 use crate::profile::{DFTProfile, MAX_POTENTIAL};
 use crate::solver::DFTSolver;
 use crate::WeightFunctionInfo;
+use feos_core::si::{
+    Density, Dimensionless, Energy, Length, MolarEnergy, MolarWeight, Quantity, Temperature,
+    Volume, _Moles, _Pressure, KELVIN, RGAS,
+};
 use feos_core::{Components, Contributions, EosResult, ReferenceSystem, State, StateBuilder};
+use feos_core::{Components, Contributions, EosResult, State, StateBuilder};
 use ndarray::prelude::*;
-use ndarray::Axis as Axis_nd;
-use ndarray::RemoveAxis;
+use ndarray::{Axis as Axis_nd, RemoveAxis};
 use num_dual::linalg::LU;
 use num_dual::DualNum;
 use quantity::{Density, Dimensionless, Energy, Length, MolarEnergy, Temperature, Volume, KELVIN};
 use std::fmt::Display;
 use std::sync::Arc;
+use typenum::Diff;
 
 const POTENTIAL_OFFSET: f64 = 2.0;
 const DEFAULT_GRID_POINTS: usize = 2048;
 
+pub type _HenryCoefficient = Diff<_Moles, _Pressure>;
+pub type HenryCoefficient<T> = Quantity<T, _HenryCoefficient>;
+
 /// Parameters required to specify a 1D pore.
 pub struct Pore1D {
     pub geometry: Geometry,
@@ -144,6 +152,35 @@ where
             * Dimensionless::new(&self.profile.bulk.molefracs))
         .sum())
     }
+
+    pub fn henry_coefficient(&self) -> EosResult<HenryCoefficient<Array1<f64>>> {
+        let mut pot = self.profile.external_potential.clone();
+        // pot.outer_iter_mut()
+        //     .zip(self.profile.dft.m().iter())
+        //     .for_each(|(mut v, &m)| v /= m);
+        let exp_pot = Dimensionless::from_reduced(pot.mapv(|v| (-v).exp()));
+        let m = self.profile.dft.m().into_owned();
+        Ok(self.profile.integrate_comp(&exp_pot)
+            / (RGAS * self.profile.temperature * Dimensionless::from_reduced(m)))
+    }
+
+    pub fn ideal_gas_enthalpy_of_adsorption(&self) -> EosResult<MolarEnergy<Array1<f64>>>
+    where
+        D::Larger: Dimension<Smaller = D>,
+    {
+        let mut pot = self.profile.external_potential.clone();
+        pot.outer_iter_mut()
+            .zip(self.profile.dft.m().iter())
+            .for_each(|(mut v, &m)| v /= m);
+        let potm1 = Dimensionless::from_reduced(self.profile.external_potential.mapv(|v| 1.0 - v));
+        let exp_pot = Dimensionless::from_reduced(pot.mapv(|v| (-v).exp()));
+        let h = self
+            .profile
+            .integrate_comp(&(potm1 * &exp_pot))
+            .convert_into(self.profile.integrate_comp(&exp_pot));
+        let m = self.profile.dft.m().into_owned();
+        Ok(h * RGAS * self.profile.temperature / Dimensionless::from_reduced(m))
+    }
 }
 
 impl PoreSpecification<Ix1> for Pore1D {

diff --git a/feos-dft/src/python/adsorption/pore.rs b/feos-dft/src/python/adsorption/pore.rs
@@ -28,6 +28,7 @@ macro_rules! impl_pore {
         pub struct PyPoreProfile1D(PoreProfile1D<$func>);
 
         impl_1d_profile!(PyPoreProfile1D, [get_r, get_z]);
+        impl_pore_profile!(PyPoreProfile1D);
 
         #[pymethods]
         impl PyPore1D {
@@ -113,36 +114,14 @@ macro_rules! impl_pore {
             }
         }
 
-        #[pymethods]
-        impl PyPoreProfile1D {
-            #[getter]
-            fn get_grand_potential(&self) -> Option<Energy> {
-                self.0.grand_potential
-            }
-
-            #[getter]
-            fn get_interfacial_tension(&self) -> Option<Energy> {
-                self.0.interfacial_tension
-            }
-
-            #[getter]
-            fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult<MolarEnergy<Array1<f64>>> {
-                Ok(self.0.partial_molar_enthalpy_of_adsorption()?.into())
-            }
-
-            #[getter]
-            fn get_enthalpy_of_adsorption(&self) -> PyResult<MolarEnergy> {
-                Ok(self.0.enthalpy_of_adsorption()?.into())
-            }
-        }
-
         #[pyclass(name = "Pore2D")]
         pub struct PyPore2D(Pore2D);
 
         #[pyclass(name = "PoreProfile2D")]
         pub struct PyPoreProfile2D(PoreProfile2D<$func>);
 
         impl_2d_profile!(PyPoreProfile2D, get_x, get_y);
+        impl_pore_profile!(PyPoreProfile2D);
 
         #[pymethods]
         impl PyPore2D {
@@ -198,38 +177,12 @@ macro_rules! impl_pore {
             }
         }
 
-        #[pymethods]
-        impl PyPoreProfile2D {
-            #[getter]
-            fn get_grand_potential(&self) -> Option<Energy> {
-                self.0.grand_potential
-            }
-
-            #[getter]
-            fn get_interfacial_tension(&self) -> Option<Energy> {
-                self.0.interfacial_tension
-            }
-
-            #[getter]
-            fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult<MolarEnergy<Array1<f64>>> {
-                Ok(self.0.partial_molar_enthalpy_of_adsorption()?)
-            }
-
-            #[getter]
-            fn get_enthalpy_of_adsorption(&self) -> PyResult<MolarEnergy> {
-                Ok(self.0.enthalpy_of_adsorption()?)
-            }
-        }
-
         /// Parameters required to specify a 3D pore.
         ///
         /// Parameters
         /// ----------
         /// system_size : [SINumber; 3]
         ///     The size of the unit cell.
-        /// angles : [Angle; 3]
-        ///     The angles of the unit cell or `None` if the unit cell
-        ///     is orthorombic
         /// n_grid : [int; 3]
         ///     The number of grid points in each direction.
         /// coordinates : numpy.ndarray[float]
@@ -238,6 +191,9 @@ macro_rules! impl_pore {
         ///     The size parameters of all interaction sites.
         /// epsilon_k_ss : numpy.ndarray[float]
         ///     The energy parameter of all interaction sites.
+        /// angles : [Angle; 3], optional
+        ///     The angles of the unit cell or `None` if the unit cell
+        ///     is orthorombic
         /// potential_cutoff: float, optional
         ///     Maximum value for the external potential.
         /// cutoff_radius: SINumber, optional
@@ -254,6 +210,7 @@ macro_rules! impl_pore {
         pub struct PyPoreProfile3D(PoreProfile3D<$func>);
 
         impl_3d_profile!(PyPoreProfile3D, get_x, get_y, get_z);
+        impl_pore_profile!(PyPoreProfile3D);
 
         #[pymethods]
         impl PyPore3D {
@@ -320,9 +277,14 @@ macro_rules! impl_pore {
                 Ok(self.0.pore_volume()?.into())
             }
         }
+    };
+}
 
+#[macro_export]
+macro_rules! impl_pore_profile {
+    ($py_profile:ty) => {
         #[pymethods]
-        impl PyPoreProfile3D {
+        impl $py_profile {
             #[getter]
             fn get_grand_potential(&self) -> Option<Energy> {
                 self.0.grand_potential
@@ -334,14 +296,26 @@ macro_rules! impl_pore {
             }
 
             #[getter]
-            fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult<MolarEnergy<Array1<f64>>> {
+            fn get_partial_molar_enthalpy_of_adsorption(
+                &self,
+            ) -> PyResult<MolarEnergy<Array1<f64>>> {
                 Ok(self.0.partial_molar_enthalpy_of_adsorption()?)
             }
 
             #[getter]
             fn get_enthalpy_of_adsorption(&self) -> PyResult<MolarEnergy> {
                 Ok(self.0.enthalpy_of_adsorption()?)
             }
+
+            #[getter]
+            fn get_henry_coefficient(&self) -> PyResult<PySIArray1> {
+                Ok(self.0.henry_coefficient()?.into())
+            }
+
+            #[getter]
+            fn get_ideal_gas_enthalpy_of_adsorption(&self) -> PyResult<PySIArray1> {
+                Ok(self.0.ideal_gas_enthalpy_of_adsorption()?.into())
+            }
         }
     };
 }