Added DFTSpecification interface to PoreProfile.initialize; started i…

…mplementation of isostere calculation
feos-org · JohannesEller · Nov 13, 2021 · Nov 13, 2021 · Nov 28, 2021 · Nov 13, 2021
commit 9ebee6e91dd9c693f23a083600415759628c81f7
diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs
@@ -1,10 +1,12 @@
 //! Adsorption profiles and isotherms.
 use super::functional::{HelmholtzEnergyFunctional, DFT};
+// use super::profile::DFTSpecifications;
 use super::solver::DFTSolver;
 use feos_core::{
     Contributions, EosError, EosResult, EosUnit, EquationOfState, StateBuilder, VLEOptions,
 };
 use ndarray::{arr1, Array1, Dimension, Ix1, Ix3};
+// use quantity::si::{KB, NAV};
 use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
 use std::rc::Rc;
 
@@ -17,6 +19,34 @@ pub use pore::{Pore1D, Pore3D, PoreProfile, PoreProfile1D, PoreProfile3D, PoreSp
 const MAX_ITER_ADSORPTION_EQUILIBRIUM: usize = 50;
 const TOL_ADSORPTION_EQUILIBRIUM: f64 = 1e-8;
 
+/// Possible inputs for the temperature grid of adsorption isostere.
+// pub enum TemperatureSpecification<U> {
+//     /// Specify the minimal and maximal pressure, and the number of points
+//     Tlim {
+//         t_min: QuantityScalar<U>,
+//         t_max: QuantityScalar<U>,
+//         points: usize,
+//     },
+//     /// Provide a custom array of pressure points.
+//     Tvec(QuantityArray1<U>),
+// }
+
+// impl<U: EosUnit> TemperatureSpecification<U>
+// where
+//     QuantityScalar<U>: std::fmt::Display,
+// {
+//     fn to_vec(&self) -> EosResult<QuantityArray1<U>> {
+//         Ok(match self {
+//             Self::Tlim {
+//                 t_min,
+//                 t_max,
+//                 points,
+//             } => QuantityArray1::linspace(*t_min, *t_max, *points)?,
+//             Self::Tvec(temperature) => temperature.clone(),
+//         })
+//     }
+// }
+
 /// Possible inputs for the pressure grid of adsorption isotherms.
 pub enum PressureSpecification<U> {
     /// Specify the minimal and maximal pressure, and the number of points
@@ -247,7 +277,10 @@ where
                 .vapor()
                 .clone();
         }
-        let external_potential = pore.initialize(&bulk, None)?.profile.external_potential;
+        let external_potential = pore
+            .initialize(&bulk, None, None)?
+            .profile
+            .external_potential;
 
         for i in 0..pressure.len() {
             let mut bulk = StateBuilder::new(functional)
@@ -261,7 +294,7 @@ where
                     .vapor()
                     .clone();
             }
-            let mut p = pore.initialize(&bulk, Some(&external_potential))?;
+            let mut p = pore.initialize(&bulk, Some(&external_potential), None)?;
             let p2 = p.clone();
             if let Some(Ok(l)) = profiles.last() {
                 p.profile.density = l.profile.density.clone();
@@ -300,8 +333,8 @@ where
             .liquid()
             .build()?;
 
-        let mut vapor = pore.initialize(&vapor_bulk, None)?.solve(None)?;
-        let mut liquid = pore.initialize(&liquid_bulk, None)?.solve(solver)?;
+        let mut vapor = pore.initialize(&vapor_bulk, None, None)?.solve(None)?;
+        let mut liquid = pore.initialize(&liquid_bulk, None, None)?.solve(solver)?;
 
         // calculate initial value for the molar gibbs energy
         let nv = vapor.profile.bulk.density
@@ -372,6 +405,47 @@ where
         ))
     }
 
+    // pub fn isostere<S: PoreSpecification<U, D, F>>(
+    //     &self,
+    //     functional: &Rc<DFT<F>>,
+    //     initial_pressure: QuantityScalar<U>,
+    //     temperature: &TemperatureSpecification<U>,
+    //     pore: &S,
+    //     molefracs: Option<&Array1<f64>>,
+    //     solver: Option<&DFTSolver>,
+    // ) -> EosResult<Adsorption<U, D, F>> {
+    //     let temperature = temperature.to_vec()?;
+    //     let moles =
+    //         functional.validate_moles(molefracs.map(|x| x * U::reference_moles()).as_ref())?;
+    //     let mut profiles: Vec<EosResult<PoreProfile<U, D, F>>> =
+    //         Vec::with_capacity(temperature.len());
+
+    //     // Create the initial bulk state at the start of the isostere
+    //     let mut initial_bulk = StateBuilder::new(functional)
+    //         .temperature(temperature.get(0))
+    //         .pressure(initial_pressure)
+    //         .moles(&moles)
+    //         .build()?;
+    //     if functional.components() > 1 && !initial_bulk.is_stable(VLEOptions::default())? {
+    //         initial_bulk = initial_bulk
+    //             .tp_flash(None, VLEOptions::default(), None)?
+    //             .vapor()
+    //             .clone();
+    //     }
+
+    //     // Initial PoreProfile for given mu,V,T
+    //     let initial_profile = pore.initialize(&initial_bulk, None, None)?.solve(solver);
+
+    //     // Specify the number of particles for the isostere calculation
+    //     let n = initial_profile?.profile.moles().to_reduced(NAV)?;
+    //     let specification = Rc::new(DFTSpecifications::Moles { moles: n });
+
+    //     for i in 1..temperature.len() {}
+    //         let dummy_state =
+
+    //     unimplemented!()
+    // }
+
     pub fn pressure(&self) -> QuantityArray1<U> {
         QuantityArray1::from_shape_fn(self.0.len(), |i| match &self.0[i] {
             Ok(p) => {

diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs
@@ -2,7 +2,7 @@ use crate::adsorption::{ExternalPotential, FluidParameters};
 use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, DFT};
 use crate::geometry::{Axis, AxisGeometry, Grid};
-use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
+use crate::profile::{DFTProfile, DFTSpecifications, CUTOFF_RADIUS, MAX_POTENTIAL};
 use crate::solver::DFTSolver;
 use feos_core::{Contributions, EosResult, EosUnit, State};
 use ndarray::prelude::*;
@@ -88,6 +88,7 @@ pub trait PoreSpecification<U, D: Dimension, F> {
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array<f64, D::Larger>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<PoreProfile<U, D, F>>;
 }
 
@@ -159,6 +160,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecificati
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array2<f64>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<PoreProfile1D<U, F>> {
         let dft = &bulk.eos;
         let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS);
@@ -195,7 +197,13 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> PoreSpecificati
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
         Ok(PoreProfile {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?,
+            profile: DFTProfile::new(
+                grid,
+                convolver,
+                bulk,
+                Some(external_potential),
+                specification,
+            )?,
             grand_potential: None,
             interfacial_tension: None,
         })
@@ -209,6 +217,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional, P: FluidParameters> PoreSpecifica
         &self,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<&Array4<f64>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<PoreProfile3D<U, F>> {
         let dft = &bulk.eos;
 
@@ -251,7 +260,13 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional, P: FluidParameters> PoreSpecifica
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
         Ok(PoreProfile {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?,
+            profile: DFTProfile::new(
+                grid,
+                convolver,
+                bulk,
+                Some(external_potential),
+                specification,
+            )?,
             grand_potential: None,
             interfacial_tension: None,
         })

diff --git a/src/interface/mod.rs b/src/interface/mod.rs
@@ -82,7 +82,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, None);
 
         Ok(Self {
-            profile: DFTProfile::new(grid, convolver, vle.vapor(), None)?,
+            profile: DFTProfile::new(grid, convolver, vle.vapor(), None, None)?,
             vle: vle.clone(),
             surface_tension: None,
             equimolar_radius: None,

diff --git a/src/profile.rs b/src/profile.rs
@@ -33,6 +33,7 @@ pub trait DFTSpecification<U, D: Dimension, F> {
 }
 
 /// Common specifications for the grand potentials in a DFT calculation.
+#[derive(Clone)]
 pub enum DFTSpecifications {
     /// DFT with specified chemical potential.
     ChemicalPotential,
@@ -49,6 +50,12 @@ pub enum DFTSpecifications {
     },
 }
 
+impl Default for DFTSpecifications {
+    fn default() -> Self {
+        Self::ChemicalPotential
+    }
+}
+
 impl DFTSpecifications {
     /// Calculate the number of particles from the profile.
     ///
@@ -198,6 +205,7 @@ where
         convolver: Rc<dyn Convolver<f64, D>>,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<Array<f64, D::Larger>>,
+        specification: Option<DFTSpecifications>,
     ) -> EosResult<Self> {
         let dft = bulk.eos.clone();
 
@@ -224,14 +232,17 @@ where
                 .assign(&(isaft.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]));
         }
 
+        // initialize DFT specification
+        let specification = specification.unwrap_or_default();
+
         Ok(Self {
             grid,
             convolver,
             dft: bulk.eos.clone(),
             temperature: bulk.temperature,
             density: density * U::reference_density(),
             chemical_potential: bulk.chemical_potential(Contributions::Total),
-            specification: Rc::new(DFTSpecifications::ChemicalPotential),
+            specification: Rc::new(specification),
             external_potential,
             bulk: bulk.clone(),
         })

diff --git a/src/python/adsorption/mod.rs b/src/python/adsorption/mod.rs
@@ -2,6 +2,7 @@ mod external_potential;
 mod pore;
 
 pub use external_potential::{PyExternalPotential, PyGeometry};
+pub use pore::PyDFTSpecification;
 
 #[macro_export]
 macro_rules! impl_adsorption {