Add ideal gas contribution to DFT calculations (#10)

prehner · web-flow · commit a7665ae5d244 · 2022-02-14T19:13:43.000+01:00
* add ideal_gas getter to HelmholtzEnergyFunctional

* Add ideal gas contribution to DFT calculations
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,8 +5,11 @@ 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
+- `HelmholtzEnergyFunctional`s can now overwrite the `ideal_gas` method to provide a non-default ideal gas contribution, that is accounted for in the calculation of the entropy, the internal energy and other properties. [#10](https://github.com/feos-org/feos-core/pull/10)
+
 ### Changed
--  Removed the `functional` field in `Pore1D` and `Pore3D`. [#9](https://github.com/feos-org/feos-core/pull/9) 
+- Removed the `functional` field in `Pore1D` and `Pore3D`. [#9](https://github.com/feos-org/feos-core/pull/9)
 
 ### Fixed
 - Fixed the units of default values for adsorption isotherms. [#8](https://github.com/feos-org/feos-core/pull/8)
diff --git a/src/functional.rs b/src/functional.rs
@@ -4,23 +4,28 @@ use crate::ideal_chain_contribution::IdealChainContribution;
 use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape};
 use feos_core::{
     Contributions, EosResult, EosUnit, EquationOfState, HelmholtzEnergy, HelmholtzEnergyDual,
-    MolarWeight, StateHD,
+    IdealGasContribution, IdealGasContributionDual, MolarWeight, StateHD,
 };
 use ndarray::*;
 use num_dual::*;
 use petgraph::graph::{Graph, UnGraph};
 use petgraph::visit::EdgeRef;
 use petgraph::Directed;
 use quantity::{QuantityArray, QuantityArray1, QuantityScalar};
+use std::fmt;
 use std::ops::{AddAssign, MulAssign};
 use std::rc::Rc;
 
 /// Wrapper struct for the [HelmholtzEnergyFunctional] trait.
 #[derive(Clone)]
 pub struct DFT<T> {
+    /// Helmholtz energy functional
     pub functional: T,
+    /// map segment -> component
     pub component_index: Array1<usize>,
+    /// chain lengths of individual components
     pub m: Array1<f64>,
+    /// ideal chain contribution
     pub ideal_chain_contribution: IdealChainContribution,
 }
 
@@ -54,6 +59,19 @@ impl<T: MolarWeight<U>, U: EosUnit> MolarWeight<U> for DFT<T> {
     }
 }
 
+struct DefaultIdealGasContribution();
+impl<D: DualNum<f64>> IdealGasContributionDual<D> for DefaultIdealGasContribution {
+    fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1<D> {
+        Array1::zeros(components)
+    }
+}
+
+impl fmt::Display for DefaultIdealGasContribution {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "Ideal gas (default)")
+    }
+}
+
 impl<T: HelmholtzEnergyFunctional> EquationOfState for DFT<T> {
     fn components(&self) -> usize {
         self.component_index[self.component_index.len() - 1] + 1
@@ -107,6 +125,10 @@ impl<T: HelmholtzEnergyFunctional> EquationOfState for DFT<T> {
         ));
         res
     }
+
+    fn ideal_gas(&self) -> &dyn IdealGasContribution {
+        self.functional.ideal_gas()
+    }
 }
 
 /// A general Helmholtz energy functional.
@@ -125,6 +147,18 @@ pub trait HelmholtzEnergyFunctional: Sized {
     /// equation of state anyways).
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64;
 
+    /// Return the ideal gas contribution.
+    ///
+    /// Per default this function returns an ideal gas contribution
+    /// in which the de Broglie wavelength is 1 for every component.
+    /// Therefore, the correct ideal gas pressure is obtained even
+    /// with no explicit ideal gas term. If a more detailed model is
+    /// required (e.g. for the calculation of internal energies) this
+    /// function has to be overwritten.
+    fn ideal_gas(&self) -> &dyn IdealGasContribution {
+        &DefaultIdealGasContribution()
+    }
+
     /// Overwrite this, if the functional consists of heterosegmented chains.
     fn bond_lengths(&self, _temperature: f64) -> UnGraph<(), f64> {
         Graph::with_capacity(0, 0)
@@ -139,6 +173,7 @@ pub trait HelmholtzEnergyFunctional: Sized {
 }
 
 impl<T: HelmholtzEnergyFunctional> DFT<T> {
+    /// Calculate the grand potential density $\omega$.
     pub fn grand_potential_density<U, D>(
         &self,
         temperature: QuantityScalar<U>,
@@ -169,12 +204,49 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         Ok(f * t * U::reference_pressure())
     }
 
+    pub(crate) fn ideal_gas_contribution<D>(
+        &self,
+        temperature: f64,
+        density: &Array<f64, D::Larger>,
+    ) -> Array<f64, D>
+    where
+        D: Dimension,
+        D::Larger: Dimension<Smaller = D>,
+    {
+        let n = self.components();
+        let ig = self.functional.ideal_gas();
+        let lambda = ig.de_broglie_wavelength(temperature, n);
+        let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
+        for (i, rhoi) in density.outer_iter().enumerate() {
+            phi += &rhoi.mapv(|rhoi| (rhoi.ln() + lambda[i] - 1.0) * rhoi);
+        }
+        phi * temperature
+    }
+
+    fn ideal_gas_contribution_dual<D>(
+        &self,
+        temperature: Dual64,
+        density: &Array<f64, D::Larger>,
+    ) -> Array<Dual64, D>
+    where
+        D: Dimension,
+        D::Larger: Dimension<Smaller = D>,
+    {
+        let n = self.components();
+        let ig = self.functional.ideal_gas();
+        let lambda = ig.de_broglie_wavelength(temperature, n);
+        let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
+        for (i, rhoi) in density.outer_iter().enumerate() {
+            phi += &rhoi.mapv(|rhoi| (lambda[i] + rhoi.ln() - 1.0) * rhoi);
+        }
+        phi * temperature
+    }
+
     fn intrinsic_helmholtz_energy_density<D, N>(
         &self,
         temperature: N,
         density: &Array<f64, D::Larger>,
         convolver: &Rc<dyn Convolver<N, D>>,
-        contributions: Contributions,
     ) -> EosResult<Array<N, D>>
     where
         N: DualNum<f64> + ScalarOperand,
@@ -187,7 +259,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         let functional_contributions = self.functional.contributions();
         let mut helmholtz_energy_density: Array<N, D> = self
             .ideal_chain_contribution
-            .calculate_helmholtz_energy_density(&density.mapv(N::from), contributions)?;
+            .calculate_helmholtz_energy_density(&density.mapv(N::from))?;
         for (c, wd) in functional_contributions.iter().zip(weighted_densities) {
             let nwd = wd.shape()[0];
             let ngrid = wd.len() / nwd;
@@ -203,6 +275,9 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         Ok(helmholtz_energy_density * temperature)
     }
 
+    /// Calculate the entropy density $s$.
+    ///
+    /// Untested with heterosegmented functionals.
     pub fn entropy_density<D>(
         &self,
         temperature: f64,
@@ -215,21 +290,26 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         D::Larger: Dimension<Smaller = D>,
     {
         let temperature_dual = Dual64::from(temperature).derive();
-        let helmholtz_energy_density = self.intrinsic_helmholtz_energy_density(
-            temperature_dual,
-            density,
-            convolver,
-            contributions,
-        )?;
+        let mut helmholtz_energy_density =
+            self.intrinsic_helmholtz_energy_density(temperature_dual, density, convolver)?;
+        match contributions {
+            Contributions::Total => {
+                helmholtz_energy_density += &self.ideal_gas_contribution_dual::<D>(temperature_dual, density);
+            },
+            Contributions::ResidualP|Contributions::IdealGas => panic!("Entropy density can only be calculated for Contributions::Residual or Contributions::Total"),
+            Contributions::Residual => (),
+        }
         Ok(helmholtz_energy_density.mapv(|f| -f.eps[0]))
     }
 
+    /// Calculate the individual contributions to the entropy density.
+    ///
+    /// Untested with heterosegmented functionals.
     pub fn entropy_density_contributions<D>(
         &self,
         temperature: f64,
         density: &Array<f64, D::Larger>,
         convolver: &Rc<dyn Convolver<Dual64, D>>,
-        contributions: Contributions,
     ) -> EosResult<Vec<Array<f64, D>>>
     where
         D: Dimension,
@@ -244,7 +324,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
             Vec::with_capacity(functional_contributions.len() + 1);
         helmholtz_energy_density.push(
             self.ideal_chain_contribution
-                .calculate_helmholtz_energy_density(&density.mapv(Dual64::from), contributions)?,
+                .calculate_helmholtz_energy_density(&density.mapv(Dual64::from))?,
         );
 
         for (c, wd) in functional_contributions.iter().zip(weighted_densities) {
@@ -265,6 +345,9 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
             .collect())
     }
 
+    /// Calculate the internal energy density $u$.
+    ///
+    /// Untested with heterosegmented functionals.
     pub fn internal_energy_density<D>(
         &self,
         temperature: f64,
@@ -278,18 +361,22 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         D::Larger: Dimension<Smaller = D>,
     {
         let temperature_dual = Dual64::from(temperature).derive();
-        let helmholtz_energy_density_dual = self.intrinsic_helmholtz_energy_density(
-            temperature_dual,
-            density,
-            convolver,
-            contributions,
-        )?;
+        let mut helmholtz_energy_density_dual =
+            self.intrinsic_helmholtz_energy_density(temperature_dual, density, convolver)?;
+        match contributions {
+                Contributions::Total => {
+                    helmholtz_energy_density_dual += &self.ideal_gas_contribution_dual::<D>(temperature_dual, density);
+                },
+                Contributions::ResidualP|Contributions::IdealGas => panic!("Internal energy density can only be calculated for Contributions::Residual or Contributions::Total"),
+                Contributions::Residual => (),
+            }
         let helmholtz_energy_density = helmholtz_energy_density_dual
             .mapv(|f| f.re - f.eps[0] * temperature)
             + (external_potential * density).sum_axis(Axis(0)) * temperature;
         Ok(helmholtz_energy_density)
     }
 
+    /// Calculate the (residual) functional derivative $\frac{\delta\mathcal{F}}{\delta\rho_i(\mathbf{r})}$.
     #[allow(clippy::type_complexity)]
     pub fn functional_derivative<D>(
         &self,
@@ -325,8 +412,8 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         ))
     }
 
-    // iSAFT correction to the functional derivative
-    pub fn isaft_integrals<D>(
+    /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
+    pub fn bond_integrals<D>(
         &self,
         temperature: f64,
         functional_derivative: &Array<f64, D::Larger>,
@@ -338,13 +425,13 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
     {
         // calculate weight functions
         let bond_lengths = self.functional.bond_lengths(temperature).into_edge_type();
-        let mut isaft_weight_functions = bond_lengths.map(
+        let mut bond_weight_functions = bond_lengths.map(
             |_, _| (),
             |_, &l| WeightFunction::new_scaled(arr1(&[l]), WeightFunctionShape::Delta),
         );
         for n in bond_lengths.node_indices() {
             for e in bond_lengths.edges(n) {
-                isaft_weight_functions.add_edge(
+                bond_weight_functions.add_edge(
                     e.target(),
                     e.source(),
                     WeightFunction::new_scaled(arr1(&[*e.weight()]), WeightFunctionShape::Delta),
@@ -354,7 +441,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
 
         let expdfdrho = functional_derivative.mapv(|x| (-x).exp());
         let mut i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
-            isaft_weight_functions.map(|_, _| (), |_, _| None);
+            bond_weight_functions.map(|_, _| (), |_, _| None);
 
         let bonds = i_graph.edge_count();
         let mut calc = 0;
@@ -384,9 +471,8 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
                                 .to_owned(),
                             |acc: Array<f64, D>, e| acc * e.weight().as_ref().unwrap(),
                         );
-                        i1 = Some(
-                            convolver.convolve(i0.clone(), &isaft_weight_functions[edge.id()]),
-                        );
+                        i1 =
+                            Some(convolver.convolve(i0.clone(), &bond_weight_functions[edge.id()]));
                         break 'nodes;
                     }
                 }
diff --git a/src/ideal_chain_contribution.rs b/src/ideal_chain_contribution.rs
@@ -1,4 +1,4 @@
-use feos_core::{Contributions, EosResult, EosUnit, HelmholtzEnergyDual, StateHD};
+use feos_core::{EosResult, EosUnit, HelmholtzEnergyDual, StateHD};
 use ndarray::*;
 use num_dual::DualNum;
 use quantity::{QuantityArray, QuantityScalar};
@@ -48,22 +48,15 @@ impl IdealChainContribution {
     pub fn calculate_helmholtz_energy_density<D, N>(
         &self,
         density: &Array<N, D::Larger>,
-        contributions: Contributions,
     ) -> EosResult<Array<N, D>>
     where
         D: Dimension,
         D::Larger: Dimension<Smaller = D>,
         N: DualNum<f64>,
     {
         let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
-        let m = match contributions {
-            Contributions::Total => self.m.clone(),
-            Contributions::Residual => self.m.clone() - 1.0,
-            Contributions::IdealGas => Array::ones(density.shape()[0]),
-            Contributions::ResidualP => unreachable!(),
-        };
         for (i, rhoi) in density.outer_iter().enumerate() {
-            phi = phi + rhoi.mapv(|rhoi| (rhoi.ln() - 1.0) * m[i] * rhoi);
+            phi += &rhoi.mapv(|rhoi| (rhoi.ln() - 1.0) * (self.m[i] - 1.0) * rhoi);
         }
         Ok(phi)
     }
@@ -72,18 +65,13 @@ impl IdealChainContribution {
         &self,
         temperature: QuantityScalar<U>,
         density: &QuantityArray<U, D::Larger>,
-        contributions: Contributions,
     ) -> EosResult<QuantityArray<U, D>>
     where
         D: Dimension,
         D::Larger: Dimension<Smaller = D>,
     {
         let rho = density.to_reduced(U::reference_density())?;
         let t = temperature.to_reduced(U::reference_temperature())?;
-        Ok(
-            self.calculate_helmholtz_energy_density(&rho, contributions)?
-                * t
-                * U::reference_pressure(),
-        )
+        Ok(self.calculate_helmholtz_energy_density(&rho)? * t * U::reference_pressure())
     }
 }
diff --git a/src/pdgt.rs b/src/pdgt.rs
@@ -169,11 +169,14 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         }
         delta_omega += &self
             .ideal_chain_contribution
-            .helmholtz_energy_density::<_, Ix1>(
-                vle.vapor().temperature,
-                &density,
-                Contributions::Total,
-            )?;
+            .helmholtz_energy_density::<_, Ix1>(vle.vapor().temperature, &density)?;
+
+        let t = vle
+            .vapor()
+            .temperature
+            .to_reduced(U::reference_temperature())?;
+        let rho = density.to_reduced(U::reference_density())?;
+        delta_omega += &(self.ideal_gas_contribution::<Ix1>(t, &rho) * U::reference_pressure());
 
         // calculate excess grand potential density
         let mu = vle.vapor().chemical_potential(Contributions::Total);
diff --git a/src/profile.rs b/src/profile.rs
