From 5f83f5feee7526aea3382e3f18f67292b2024a30 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 14 Feb 2022 15:44:12 +0100 Subject: [PATCH 1/2] add ideal_gas getter to HelmholtzEnergyFunctional --- src/functional.rs | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/functional.rs b/src/functional.rs index acd3033..580b05e 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -4,7 +4,7 @@ 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::*; @@ -12,6 +12,7 @@ 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; @@ -54,6 +55,19 @@ impl, U: EosUnit> MolarWeight for DFT { } } +struct DefaultIdealGasContribution(); +impl> IdealGasContributionDual for DefaultIdealGasContribution { + fn de_broglie_wavelength(&self, _: D, components: usize) -> Array1 { + Array1::zeros(components) + } +} + +impl fmt::Display for DefaultIdealGasContribution { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Ideal gas (default)") + } +} + impl EquationOfState for DFT { fn components(&self) -> usize { self.component_index[self.component_index.len() - 1] + 1 @@ -107,6 +121,10 @@ impl EquationOfState for DFT { )); res } + + fn ideal_gas(&self) -> &dyn IdealGasContribution { + self.functional.ideal_gas() + } } /// A general Helmholtz energy functional. @@ -125,6 +143,18 @@ pub trait HelmholtzEnergyFunctional: Sized { /// equation of state anyways). fn compute_max_density(&self, moles: &Array1) -> 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) From 22b1e2ca12db8079f188ace801bf68379512108d Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 14 Feb 2022 18:35:59 +0100 Subject: [PATCH 2/2] Add ideal gas contribution to DFT calculations --- CHANGELOG.md | 5 +- src/functional.rs | 104 ++++++++++++++++++++++++-------- src/ideal_chain_contribution.rs | 18 +----- src/pdgt.rs | 13 ++-- src/profile.rs | 41 +++++++++---- 5 files changed, 125 insertions(+), 56 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b58ab1b..c14dd23 100644 --- 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 index 580b05e..4578f96 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -19,9 +19,13 @@ use std::rc::Rc; /// Wrapper struct for the [HelmholtzEnergyFunctional] trait. #[derive(Clone)] pub struct DFT { + /// Helmholtz energy functional pub functional: T, + /// map segment -> component pub component_index: Array1, + /// chain lengths of individual components pub m: Array1, + /// ideal chain contribution pub ideal_chain_contribution: IdealChainContribution, } @@ -169,6 +173,7 @@ pub trait HelmholtzEnergyFunctional: Sized { } impl DFT { + /// Calculate the grand potential density $\omega$. pub fn grand_potential_density( &self, temperature: QuantityScalar, @@ -199,12 +204,49 @@ impl DFT { Ok(f * t * U::reference_pressure()) } + pub(crate) fn ideal_gas_contribution( + &self, + temperature: f64, + density: &Array, + ) -> Array + where + D: Dimension, + D::Larger: Dimension, + { + 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( + &self, + temperature: Dual64, + density: &Array, + ) -> Array + where + D: Dimension, + D::Larger: Dimension, + { + 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( &self, temperature: N, density: &Array, convolver: &Rc>, - contributions: Contributions, ) -> EosResult> where N: DualNum + ScalarOperand, @@ -217,7 +259,7 @@ impl DFT { let functional_contributions = self.functional.contributions(); let mut helmholtz_energy_density: Array = 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; @@ -233,6 +275,9 @@ impl DFT { Ok(helmholtz_energy_density * temperature) } + /// Calculate the entropy density $s$. + /// + /// Untested with heterosegmented functionals. pub fn entropy_density( &self, temperature: f64, @@ -245,21 +290,26 @@ impl DFT { D::Larger: Dimension, { 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::(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( &self, temperature: f64, density: &Array, convolver: &Rc>, - contributions: Contributions, ) -> EosResult>> where D: Dimension, @@ -274,7 +324,7 @@ impl DFT { 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) { @@ -295,6 +345,9 @@ impl DFT { .collect()) } + /// Calculate the internal energy density $u$. + /// + /// Untested with heterosegmented functionals. pub fn internal_energy_density( &self, temperature: f64, @@ -308,18 +361,22 @@ impl DFT { D::Larger: Dimension, { 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::(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( &self, @@ -355,8 +412,8 @@ impl DFT { )) } - // iSAFT correction to the functional derivative - pub fn isaft_integrals( + /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$ + pub fn bond_integrals( &self, temperature: f64, functional_derivative: &Array, @@ -368,13 +425,13 @@ impl DFT { { // 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), @@ -384,7 +441,7 @@ impl DFT { let expdfdrho = functional_derivative.mapv(|x| (-x).exp()); let mut i_graph: Graph<_, Option>, Directed> = - isaft_weight_functions.map(|_, _| (), |_, _| None); + bond_weight_functions.map(|_, _| (), |_, _| None); let bonds = i_graph.edge_count(); let mut calc = 0; @@ -414,9 +471,8 @@ impl DFT { .to_owned(), |acc: Array, 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 index 67d04d1..7fa9a00 100644 --- 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,7 +48,6 @@ impl IdealChainContribution { pub fn calculate_helmholtz_energy_density( &self, density: &Array, - contributions: Contributions, ) -> EosResult> where D: Dimension, @@ -56,14 +55,8 @@ impl IdealChainContribution { N: DualNum, { 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,7 +65,6 @@ impl IdealChainContribution { &self, temperature: QuantityScalar, density: &QuantityArray, - contributions: Contributions, ) -> EosResult> where D: Dimension, @@ -80,10 +72,6 @@ impl IdealChainContribution { { 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 index f11add5..30bd262 100644 --- a/src/pdgt.rs +++ b/src/pdgt.rs @@ -169,11 +169,14 @@ impl DFT { } 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::(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 index 12a400c..e779ac0 100644 --- a/src/profile.rs +++ b/src/profile.rs @@ -212,8 +212,8 @@ where // intitialize density let t = bulk.temperature.to_reduced(U::reference_temperature())?; - let isaft = dft - .isaft_integrals(t, &external_potential, &convolver) + 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()); @@ -221,7 +221,7 @@ where for (s, &c) in dft.component_index.iter().enumerate() { density .index_axis_mut(Axis_nd(0), s) - .assign(&(isaft.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c])); + .assign(&(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c])); } Ok(Self { @@ -348,10 +348,16 @@ where // Read from profile let temperature = self.temperature.to_reduced(U::reference_temperature())?; let density = self.density.to_reduced(U::reference_density())?; + let lambda_de_broglie = self + .dft + .functional + .ideal_gas() + .de_broglie_wavelength(temperature, self.bulk.eos.components()); let mu_comp = self .chemical_potential .to_reduced(U::reference_molar_energy())? - / temperature; + / temperature + - lambda_de_broglie; let chemical_potential = self.dft.component_index.mapv(|i| mu_comp[i]); let mut bulk = self.bulk.clone(); @@ -383,11 +389,18 @@ where log: bool, ) -> EosResult<()> { // Update bulk state + let lambda_de_broglie = self + .dft + .functional + .ideal_gas() + .de_broglie_wavelength(temperature, bulk.eos.components()); let mut mu_comp = Array::zeros(bulk.eos.components()); for (s, &c) in self.dft.component_index.iter().enumerate() { mu_comp[c] = chemical_potential[s]; } - bulk.update_chemical_potential(&(mu_comp * temperature * U::reference_molar_energy()))?; + bulk.update_chemical_potential( + &((mu_comp + lambda_de_broglie) * temperature * U::reference_molar_energy()), + )?; // calculate intrinsic functional derivative let (_, mut dfdrho) = @@ -397,10 +410,10 @@ where // calculate total functional derivative dfdrho += &self.external_potential; - // calculate isaft integrals - let isaft = self + // calculate bond integrals + let bonds = self .dft - .isaft_integrals(temperature, &dfdrho, &self.convolver); + .bond_integrals(temperature, &dfdrho, &self.convolver); // Euler-Lagrange equation let m = &self.dft.m; @@ -410,7 +423,7 @@ where .zip(chemical_potential.iter()) .zip(m.iter()) .zip(density.outer_iter()) - .zip(isaft.outer_iter()) + .zip(bonds.outer_iter()) .for_each(|(((((mut res, df), &mu), &m), rho), is)| { res.assign( &(if log { @@ -435,7 +448,7 @@ where let z: Array1<_> = dfdrho .outer_iter() .zip(m.iter()) - .zip(isaft.outer_iter()) + .zip(bonds.outer_iter()) .map(|((df, &m), is)| self.integrate_reduced((-&df / m).mapv(f64::exp) * is)) .collect(); let mu_spec = @@ -460,10 +473,16 @@ where // Read from profile let temperature = self.temperature.to_reduced(U::reference_temperature())?; let mut density = self.density.to_reduced(U::reference_density())?; + let lambda_de_broglie = self + .dft + .functional + .ideal_gas() + .de_broglie_wavelength(temperature, self.bulk.eos.components()); let mut mu_comp = self .chemical_potential .to_reduced(U::reference_molar_energy())? - / temperature; + / temperature + - lambda_de_broglie; let mut chemical_potential = self.dft.component_index.mapv(|i| mu_comp[i]); let mut bulk = self.bulk.clone();