Add density as additional parameter to DFTProfile::new() (#24)

* Add density as additional parameter to DFTProfile::new() * Fix calculation of pore volume * update changelog
feos-org · Apr 7, 2022 · 30300da · 30300da
1 parent 82a6a16
commit 30300da
Show file tree

Hide file tree

Showing 8 changed files with 57 additions and 33 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,12 +6,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
-## [0.2.0] - 2022-03-??
+## [0.2.0] - 2022-04-??
 ### Changed
 - Renamed `AxisGeometry` to `Geometry`. [#19](https://github.com/feos-org/feos-dft/pull/19)
 - Removed `PyGeometry` and `PyFMTVersion` in favor of a simpler implementation using `PyO3`'s new `#[pyclass]` for fieldless enums feature. [#19](https://github.com/feos-org/feos-dft/pull/19)
 - `DFTSolver` now uses `Verbosity` instead of a `bool` to control its output. [#19](https://github.com/feos-org/feos-dft/pull/19)
 - `SurfaceTensionDiagram` now uses the new `StateVec` struct to access properties of the bulk phases. [#19](https://github.com/feos-org/feos-dft/pull/19)
+- `Pore1D::initialize` and `Pore3D::initialize` now accept initial values for the density profiles as optional arguments. [#24](https://github.com/feos-org/feos-dft/pull/24)
 
 ### Packaging
 - Updated `pyo3` and `numpy` dependencies to 0.16.

diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs
@@ -265,7 +265,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)
@@ -279,11 +282,14 @@ where
                     .vapor()
                     .clone();
             }
-            let mut p = pore.initialize(&bulk, Some(&external_potential))?;
-            let p2 = p.clone();
-            if let Some(Ok(l)) = profiles.last() {
-                p.profile.density = l.profile.density.clone();
-            }
+            let old_density = if let Some(Ok(l)) = profiles.last() {
+                Some(&l.profile.density)
+            } else {
+                None
+            };
+
+            let p = pore.initialize(&bulk, old_density, Some(&external_potential))?;
+            let p2 = pore.initialize(&bulk, None, Some(&external_potential))?;
             profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));
         }
 
@@ -318,8 +324,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

diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs
@@ -10,7 +10,7 @@ use ndarray::prelude::*;
 use ndarray::Axis as Axis_nd;
 use ndarray::Zip;
 use ndarray_stats::QuantileExt;
-use quantity::{QuantityArray2, QuantityScalar};
+use quantity::{QuantityArray, QuantityArray2, QuantityArray4, QuantityScalar};
 use std::rc::Rc;
 
 const POTENTIAL_OFFSET: f64 = 2.0;
@@ -83,6 +83,7 @@ pub trait PoreSpecification<U: EosUnit, D: Dimension> {
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<U, DFT<F>>,
+        density: Option<&QuantityArray<U, D::Larger>>,
         external_potential: Option<&Array<f64, D::Larger>>,
     ) -> EosResult<PoreProfile<U, D, F>>;
 
@@ -96,9 +97,9 @@ pub trait PoreSpecification<U: EosUnit, D: Dimension> {
     {
         let bulk = StateBuilder::new(&Rc::new(Helium::new()))
             .temperature(298.0 * U::reference_temperature())
-            .volume(U::reference_volume())
+            .density(U::reference_density())
             .build()?;
-        let pore = self.initialize(&bulk, None)?;
+        let pore = self.initialize(&bulk, None, None)?;
         let pot = pore
             .profile
             .external_potential
@@ -174,6 +175,7 @@ impl<U: EosUnit> PoreSpecification<U, Ix1> for Pore1D<U> {
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<U, DFT<F>>,
+        density: Option<&QuantityArray2<U>>,
         external_potential: Option<&Array2<f64>>,
     ) -> EosResult<PoreProfile1D<U, F>> {
         let dft = &bulk.eos;
@@ -211,7 +213,7 @@ impl<U: EosUnit> PoreSpecification<U, Ix1> for Pore1D<U> {
         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), density)?,
             grand_potential: None,
             interfacial_tension: None,
         })
@@ -226,6 +228,7 @@ impl<U: EosUnit> PoreSpecification<U, Ix3> for Pore3D<U> {
     fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
         &self,
         bulk: &State<U, DFT<F>>,
+        density: Option<&QuantityArray4<U>>,
         external_potential: Option<&Array4<f64>>,
     ) -> EosResult<PoreProfile3D<U, F>> {
         let dft = &bulk.eos;
@@ -269,7 +272,7 @@ impl<U: EosUnit> PoreSpecification<U, Ix3> for Pore3D<U> {
         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), density)?,
             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
@@ -189,6 +189,7 @@ where
         convolver: Rc<dyn Convolver<f64, D>>,
         bulk: &State<U, DFT<F>>,
         external_potential: Option<Array<f64, D::Larger>>,
+        density: Option<&QuantityArray<U, D::Larger>>,
     ) -> EosResult<Self> {
         let dft = bulk.eos.clone();
 
@@ -201,26 +202,31 @@ where
             Array::zeros(n_grid).into_dimensionality().unwrap()
         });
 
-        // intitialize density
-        let t = bulk.temperature.to_reduced(U::reference_temperature())?;
-        let bonds = dft
-            .bond_integrals(t, &external_potential, &convolver)
-            .mapv(f64::abs)
-            * (-&external_potential).mapv(f64::exp);
-        let mut density = Array::zeros(external_potential.raw_dim());
-        let bulk_density = bulk.partial_density.to_reduced(U::reference_density())?;
-        for (s, &c) in dft.component_index.iter().enumerate() {
-            density
-                .index_axis_mut(Axis_nd(0), s)
-                .assign(&(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]));
-        }
+        // initialize density
+        let density = if let Some(density) = density {
+            density.clone()
+        } else {
+            let t = bulk.temperature.to_reduced(U::reference_temperature())?;
+            let bonds = dft
+                .bond_integrals(t, &external_potential, &convolver)
+                .mapv(f64::abs)
+                * (-&external_potential).mapv(f64::exp);
+            let mut density = Array::zeros(external_potential.raw_dim());
+            let bulk_density = bulk.partial_density.to_reduced(U::reference_density())?;
+            for (s, &c) in dft.component_index.iter().enumerate() {
+                density.index_axis_mut(Axis_nd(0), s).assign(
+                    &(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]),
+                );
+            }
+            density * U::reference_density()
+        };
 
         Ok(Self {
             grid,
             convolver,
             dft: bulk.eos.clone(),
             temperature: bulk.temperature,
-            density: density * U::reference_density(),
+            density,
             specification: Rc::new(DFTSpecifications::ChemicalPotential),
             external_potential,
             bulk: bulk.clone(),

diff --git a/src/python/adsorption/pore.rs b/src/python/adsorption/pore.rs
@@ -55,6 +55,8 @@ macro_rules! impl_pore {
             /// ----------
             /// bulk : State
             ///     The bulk state in equilibrium with the pore.
+            /// density : SIArray2, optional
+            ///     Initial values for the density profile.
             /// external_potential : numpy.ndarray[float], optional
             ///     The external potential in the pore. Used to
             ///     save computation time in the case of costly
@@ -63,14 +65,16 @@ macro_rules! impl_pore {
             /// Returns
             /// -------
             /// PoreProfile1D
-            #[pyo3(text_signature = "($self, bulk, external_potential=None)")]
+            #[pyo3(text_signature = "($self, bulk, density=None, external_potential=None)")]
             fn initialize(
                 &self,
                 bulk: &PyState,
+                density: Option<PySIArray2>,
                 external_potential: Option<&PyArray2<f64>>,
             ) -> PyResult<PyPoreProfile1D> {
                 Ok(PyPoreProfile1D(self.0.initialize(
                     &bulk.0,
+                    density.as_deref(),
                     external_potential.map(|e| e.to_owned_array()).as_ref(),
                 )?))
             }
@@ -156,6 +160,8 @@ macro_rules! impl_pore {
             /// ----------
             /// bulk : State
             ///     The bulk state in equilibrium with the pore.
+            /// density : SIArray4, optional
+            ///     Initial values for the density profile.
             /// external_potential : numpy.ndarray[float], optional
             ///     The external potential in the pore. Used to
             ///     save computation time in the case of costly
@@ -164,14 +170,16 @@ macro_rules! impl_pore {
             /// Returns
             /// -------
             /// PoreProfile3D
-            #[pyo3(text_signature = "($self, bulk, external_potential=None)")]
+            #[pyo3(text_signature = "($self, bulk, density=None, external_potential=None)")]
             fn initialize(
                 &self,
                 bulk: &PyState,
+                density: Option<PySIArray4>,
                 external_potential: Option<&PyArray4<f64>>,
             ) -> PyResult<PyPoreProfile3D> {
                 Ok(PyPoreProfile3D(self.0.initialize(
                     &bulk.0,
+                    density.as_deref(),
                     external_potential.map(|e| e.to_owned_array()).as_ref(),
                 )?))
             }

diff --git a/src/solvation/mod.rs b/src/solvation/mod.rs
@@ -120,7 +120,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> SolvationProfil
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
         Ok(Self {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?,
+            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None)?,
             grand_potential: None,
             solvation_free_energy: None,
         })

diff --git a/src/solvation/pair_correlation.rs b/src/solvation/pair_correlation.rs
@@ -59,7 +59,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<U
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
         Ok(Self {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?,
+            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None)?,
             pair_correlation_function: None,
             self_solvation_free_energy: None,
             structure_factor: None,