From ad478548e86a94ee8cdf5ab22c6f7b3dba60afb1 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sun, 3 Jul 2022 15:02:02 +0200 Subject: [PATCH 1/9] New workspace crate for reusable implementations of SAFT --- Cargo.toml | 12 +- feos-dft/src/lib.rs | 1 - feos-saft/Cargo.toml | 29 ++++ feos-saft/docs-header.html | 15 ++ feos-saft/license-apache | 1 + feos-saft/license-mit | 1 + .../src/fundamental_measure_theory.rs | 131 ++++++--------- feos-saft/src/hard_sphere.rs | 153 ++++++++++++++++++ feos-saft/src/lib.rs | 7 + src/gc_pcsaft/dft/hard_chain.rs | 3 +- src/gc_pcsaft/dft/mod.rs | 12 +- src/gc_pcsaft/eos/association.rs | 10 +- src/gc_pcsaft/eos/dispersion.rs | 5 +- src/gc_pcsaft/eos/hard_chain.rs | 6 +- src/gc_pcsaft/eos/mod.rs | 3 +- src/gc_pcsaft/eos/parameter.rs | 16 ++ src/gc_pcsaft/eos/polar.rs | 1 + src/pcsaft/dft/association.rs | 1 + src/pcsaft/dft/dispersion.rs | 1 + src/pcsaft/dft/hard_chain.rs | 4 +- src/pcsaft/dft/mod.rs | 21 +-- src/pcsaft/dft/polar.rs | 1 + src/pcsaft/dft/pure_saft_functional.rs | 2 +- src/pcsaft/eos/association.rs | 1 + src/pcsaft/eos/dispersion.rs | 1 + src/pcsaft/eos/hard_chain.rs | 11 +- src/pcsaft/eos/mod.rs | 3 +- src/pcsaft/eos/polar.rs | 1 + src/pcsaft/parameters.rs | 15 ++ src/python/dft.rs | 2 +- tests/pcsaft/dft.rs | 2 +- 31 files changed, 332 insertions(+), 140 deletions(-) create mode 100644 feos-saft/Cargo.toml create mode 100644 feos-saft/docs-header.html create mode 120000 feos-saft/license-apache create mode 120000 feos-saft/license-mit rename {feos-dft => feos-saft}/src/fundamental_measure_theory.rs (75%) create mode 100644 feos-saft/src/hard_sphere.rs create mode 100644 feos-saft/src/lib.rs diff --git a/Cargo.toml b/Cargo.toml index effc378ed..e7829b771 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -16,7 +16,7 @@ categories = ["science"] features = ["all_models"] [workspace] -members = ["feos-core", "feos-dft"] +members = ["feos-core", "feos-dft", "feos-saft"] [lib] crate-type = ["rlib", "cdylib"] @@ -26,7 +26,7 @@ quantity = "0.5" num-dual = "0.5" feos-core = { version = "0.3", path = "feos-core" } feos-dft = { version = "0.2", path = "feos-dft", optional = true } -#feos-pets = { version = "0.1", features = ["python"] } +feos-saft = { version = "0.1", path = "feos-saft", optional = true } numpy = { version = "0.16", optional = true } ndarray = { version = "0.15", features = ["approx"] } petgraph = { version = "0.6", optional = true } @@ -48,11 +48,11 @@ approx = "0.4" [features] default = [] -dft = ["feos-dft", "petgraph"] +dft = ["feos-dft", "petgraph", "feos-saft?/dft"] estimator = [] -pcsaft = [] -gc_pcsaft = [] +pcsaft = ["feos-saft"] +gc_pcsaft = ["feos-saft"] uvtheory = ["lazy_static"] pets = [] -python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python"] +python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python", "feos-saft?/python"] all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets"] diff --git a/feos-dft/src/lib.rs b/feos-dft/src/lib.rs index c6020cc70..1addd3bd7 100644 --- a/feos-dft/src/lib.rs +++ b/feos-dft/src/lib.rs @@ -7,7 +7,6 @@ pub mod adsorption; mod convolver; mod functional; mod functional_contribution; -pub mod fundamental_measure_theory; mod geometry; mod ideal_chain_contribution; pub mod interface; diff --git a/feos-saft/Cargo.toml b/feos-saft/Cargo.toml new file mode 100644 index 000000000..549e80c45 --- /dev/null +++ b/feos-saft/Cargo.toml @@ -0,0 +1,29 @@ +[package] +name = "feos-saft" +version = "0.1.0" +authors = ["Gernot Bauer ", + "Philipp Rehner + + + diff --git a/feos-saft/license-apache b/feos-saft/license-apache new file mode 120000 index 000000000..67e3ce635 --- /dev/null +++ b/feos-saft/license-apache @@ -0,0 +1 @@ +../license-apache \ No newline at end of file diff --git a/feos-saft/license-mit b/feos-saft/license-mit new file mode 120000 index 000000000..137b4b11d --- /dev/null +++ b/feos-saft/license-mit @@ -0,0 +1 @@ +../license-mit \ No newline at end of file diff --git a/feos-dft/src/fundamental_measure_theory.rs b/feos-saft/src/fundamental_measure_theory.rs similarity index 75% rename from feos-dft/src/fundamental_measure_theory.rs rename to feos-saft/src/fundamental_measure_theory.rs index a721102fa..6e346b810 100644 --- a/feos-dft/src/fundamental_measure_theory.rs +++ b/feos-saft/src/fundamental_measure_theory.rs @@ -1,49 +1,22 @@ //! Helmholtz energy functionals from fundamental measure theory. -use crate::adsorption::FluidParameters; -use crate::functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use crate::functional_contribution::*; -use crate::solvation::PairPotential; -use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use feos_core::EosResult; +use feos_dft::adsorption::FluidParameters; +use feos_dft::solvation::PairPotential; +use feos_dft::{ + FunctionalContribution, FunctionalContributionDual, HelmholtzEnergyFunctional, MoleculeShape, + WeightFunction, WeightFunctionInfo, WeightFunctionShape, DFT, +}; use ndarray::*; use num_dual::DualNum; use std::f64::consts::PI; use std::fmt; use std::rc::Rc; +use crate::{HardSphereProperties, MonomerShape}; + const PI36M1: f64 = 1.0 / (36.0 * PI); const N3_CUTOFF: f64 = 1e-5; -/// Different monomer shapes for FMT. -pub enum MonomerShape { - /// For spherical monomers, the number of components. - Spherical(usize), - /// For non-spherical molecules in a homosegmented approach, the - /// chain length parameter $m$. - NonSpherical(Array1), - /// For non-spherical molecules in a heterosegmented approach, - /// the geometry factors for every segment. - Heterosegmented([Array1; 4]), -} - -/// Properties of (generalized) hard sphere systems. -pub trait FMTProperties { - fn component_index(&self) -> Array1; - fn monomer_shape>(&self, temperature: N) -> MonomerShape; - fn hs_diameter>(&self, temperature: N) -> Array1; - - fn geometry_coefficients>(&self, temperature: N) -> [Array1; 4] { - match self.monomer_shape(temperature) { - MonomerShape::Spherical(n) => { - let m = Array1::ones(n); - [m.clone(), m.clone(), m.clone(), m] - } - MonomerShape::NonSpherical(m) => [m.clone(), m.clone(), m.clone(), m], - MonomerShape::Heterosegmented(g) => g, - } - } -} - /// Different versions of fundamental measure theory #[derive(Clone, Copy)] #[cfg_attr(feature = "python", pyo3::pyclass)] @@ -80,31 +53,34 @@ impl

FMTContribution

{ } } -impl> FunctionalContributionDual for FMTContribution

{ +impl> FunctionalContributionDual + for FMTContribution

+{ fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let r = self.properties.hs_diameter(temperature) * 0.5; let [c0, c1, c2, c3] = self.properties.geometry_coefficients(temperature); match (self.version, r.len()) { (FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear, 1) => { - WeightFunctionInfo::new(self.properties.component_index(), false).extend( - vec![ - WeightFunctionShape::Delta, - WeightFunctionShape::Theta, - WeightFunctionShape::DeltaVec, - ] - .into_iter() - .zip([c2, c3.clone(), c3]) - .map(|(s, c)| WeightFunction { - prefactor: c, - kernel_radius: r.clone(), - shape: s, - }) - .collect(), - false, - ) + WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) + .extend( + vec![ + WeightFunctionShape::Delta, + WeightFunctionShape::Theta, + WeightFunctionShape::DeltaVec, + ] + .into_iter() + .zip([c2, c3.clone(), c3]) + .map(|(s, c)| WeightFunction { + prefactor: c, + kernel_radius: r.clone(), + shape: s, + }) + .collect(), + false, + ) } (FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear, _) => { - WeightFunctionInfo::new(self.properties.component_index(), false) + WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) .add( WeightFunction { prefactor: Zip::from(&c0) @@ -161,23 +137,24 @@ impl> FunctionalContributionDual for FMTCon ) } (FMTVersion::KierlikRosinberg, _) => { - WeightFunctionInfo::new(self.properties.component_index(), false).extend( - vec![ - WeightFunctionShape::KR0, - WeightFunctionShape::KR1, - WeightFunctionShape::Delta, - WeightFunctionShape::Theta, - ] - .into_iter() - .zip(self.properties.geometry_coefficients(temperature)) - .map(|(s, c)| WeightFunction { - prefactor: c, - kernel_radius: r.clone(), - shape: s, - }) - .collect(), - true, - ) + WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) + .extend( + vec![ + WeightFunctionShape::KR0, + WeightFunctionShape::KR1, + WeightFunctionShape::Delta, + WeightFunctionShape::Theta, + ] + .into_iter() + .zip(self.properties.geometry_coefficients(temperature)) + .map(|(s, c)| WeightFunction { + prefactor: c, + kernel_radius: r.clone(), + shape: s, + }) + .collect(), + true, + ) } } } @@ -275,7 +252,7 @@ impl> FunctionalContributionDual for FMTCon } } -impl fmt::Display for FMTContribution

{ +impl fmt::Display for FMTContribution

{ fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let ver = match self.version { FMTVersion::WhiteBear => "WB", @@ -286,15 +263,11 @@ impl fmt::Display for FMTContribution

{ } } -struct HardSphereProperties { +struct HardSphereParameters { sigma: Array1, } -impl FMTProperties for HardSphereProperties { - fn component_index(&self) -> Array1 { - Array1::from_shape_fn(self.sigma.len(), |i| i) - } - +impl HardSphereProperties for HardSphereParameters { fn monomer_shape(&self, _: N) -> MonomerShape { MonomerShape::Spherical(self.sigma.len()) } @@ -306,14 +279,14 @@ impl FMTProperties for HardSphereProperties { /// [HelmholtzEnergyFunctional] for hard sphere systems. pub struct FMTFunctional { - properties: Rc, + properties: Rc, contributions: Vec>, version: FMTVersion, } impl FMTFunctional { pub fn new(sigma: &Array1, version: FMTVersion) -> DFT { - let properties = Rc::new(HardSphereProperties { + let properties = Rc::new(HardSphereParameters { sigma: sigma.clone(), }); let contributions: Vec> = diff --git a/feos-saft/src/hard_sphere.rs b/feos-saft/src/hard_sphere.rs new file mode 100644 index 000000000..5bf17e16c --- /dev/null +++ b/feos-saft/src/hard_sphere.rs @@ -0,0 +1,153 @@ +use feos_core::{HelmholtzEnergyDual, StateHD}; +use ndarray::*; +use num_dual::DualNum; +use std::borrow::Cow; +use std::f64::consts::FRAC_PI_6; +use std::fmt; +use std::rc::Rc; + +/// Different monomer shapes for FMT. +pub enum MonomerShape<'a, D> { + /// For spherical monomers, the number of components. + Spherical(usize), + /// For non-spherical molecules in a homosegmented approach, the + /// chain length parameter $m$. + NonSpherical(Array1), + /// For non-spherical molecules in a heterosegmented approach, + /// the geometry factors for every segment and the component + /// index for every segment. + Heterosegmented([Array1; 4], &'a Array1), +} + +/// Properties of (generalized) hard sphere systems. +pub trait HardSphereProperties { + fn monomer_shape>(&self, temperature: D) -> MonomerShape; + fn hs_diameter>(&self, temperature: D) -> Array1; + + fn component_index(&self) -> Cow> { + match self.monomer_shape(1.0) { + MonomerShape::Spherical(n) => Cow::Owned(Array1::from_shape_fn(n, |i| i)), + MonomerShape::NonSpherical(m) => Cow::Owned(Array1::from_shape_fn(m.len(), |i| i)), + MonomerShape::Heterosegmented(_, component_index) => Cow::Borrowed(component_index), + } + } + + fn geometry_coefficients>(&self, temperature: D) -> [Array1; 4] { + match self.monomer_shape(temperature) { + MonomerShape::Spherical(n) => { + let m = Array1::ones(n); + [m.clone(), m.clone(), m.clone(), m] + } + MonomerShape::NonSpherical(m) => [m.clone(), m.clone(), m.clone(), m], + MonomerShape::Heterosegmented(g, _) => g, + } + } + + fn zeta, const N: usize>( + &self, + temperature: D, + partial_density: &Array1, + k: [i32; N], + ) -> [D; N] { + let component_index = self.component_index(); + let geometry_coefficients = self.geometry_coefficients(temperature); + let diameter = self.hs_diameter(temperature); + let mut zeta = [D::zero(); N]; + for i in 0..diameter.len() { + for (z, &k) in zeta.iter_mut().zip(k.iter()) { + *z += partial_density[component_index[i]] + * diameter[i].powi(k) + * (geometry_coefficients[k as usize][i] * FRAC_PI_6); + } + } + + zeta + } + + fn zeta_23>(&self, temperature: D, molefracs: &Array1) -> D { + let component_index = self.component_index(); + let geometry_coefficients = self.geometry_coefficients(temperature); + let diameter = self.hs_diameter(temperature); + let mut zeta: [D; 2] = [D::zero(); 2]; + for i in 0..diameter.len() { + for (k, z) in zeta.iter_mut().enumerate() { + *z += molefracs[component_index[i]] + * diameter[i].powi((k + 2) as i32) + * (geometry_coefficients[k + 2][i] * FRAC_PI_6); + } + } + + zeta[0] / zeta[1] + } +} + +pub struct HardSphere

{ + pub parameters: Rc

, +} + +impl, P: HardSphereProperties> HelmholtzEnergyDual for HardSphere

{ + fn helmholtz_energy(&self, state: &StateHD) -> D { + let p = &self.parameters; + let zeta = p.zeta(state.temperature, &state.partial_density, [0, 1, 2, 3]); + let frac_1mz3 = -(zeta[3] - 1.0).recip(); + let zeta_23 = p.zeta_23(state.temperature, &state.molefracs); + state.volume * 6.0 / std::f64::consts::PI + * (zeta[1] * zeta[2] * frac_1mz3 * 3.0 + + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23 + + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p()) + } +} + +impl

fmt::Display for HardSphere

{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Hard Sphere") + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// use crate::pcsaft::parameters::utils::{ +// butane_parameters, propane_butane_parameters, propane_parameters, +// }; +// use approx::assert_relative_eq; +// use ndarray::arr1; + +// #[test] +// fn helmholtz_energy() { +// let hs = HardSphere { +// parameters: propane_parameters(), +// }; +// let t = 250.0; +// let v = 1000.0; +// let n = 1.0; +// let s = StateHD::new(t, v, arr1(&[n])); +// let a_rust = hs.helmholtz_energy(&s); +// assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10); +// } + +// #[test] +// fn mix() { +// let c1 = HardSphere { +// parameters: propane_parameters(), +// }; +// let c2 = HardSphere { +// parameters: butane_parameters(), +// }; +// let c12 = HardSphere { +// parameters: propane_butane_parameters(), +// }; +// let t = 250.0; +// let v = 2.5e28; +// let n = 1.0; +// let s = StateHD::new(t, v, arr1(&[n])); +// let a1 = c1.helmholtz_energy(&s); +// let a2 = c2.helmholtz_energy(&s); +// let s1m = StateHD::new(t, v, arr1(&[n, 0.0])); +// let a1m = c12.helmholtz_energy(&s1m); +// let s2m = StateHD::new(t, v, arr1(&[0.0, n])); +// let a2m = c12.helmholtz_energy(&s2m); +// assert_relative_eq!(a1, a1m, epsilon = 1e-14); +// assert_relative_eq!(a2, a2m, epsilon = 1e-14); +// } +// } diff --git a/feos-saft/src/lib.rs b/feos-saft/src/lib.rs new file mode 100644 index 000000000..ce09bc930 --- /dev/null +++ b/feos-saft/src/lib.rs @@ -0,0 +1,7 @@ +mod hard_sphere; +pub use hard_sphere::{HardSphere, HardSphereProperties, MonomerShape}; + +#[cfg(feature = "dft")] +mod fundamental_measure_theory; +#[cfg(feature = "dft")] +pub use fundamental_measure_theory::{FMTContribution, FMTFunctional, FMTVersion}; diff --git a/src/gc_pcsaft/dft/hard_chain.rs b/src/gc_pcsaft/dft/hard_chain.rs index 0ff45918d..0b1628489 100644 --- a/src/gc_pcsaft/dft/hard_chain.rs +++ b/src/gc_pcsaft/dft/hard_chain.rs @@ -1,6 +1,5 @@ use super::GcPcSaftFunctionalParameters; use feos_core::EosError; -use feos_dft::fundamental_measure_theory::FMTProperties; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; @@ -27,7 +26,7 @@ impl + ScalarOperand> FunctionalContributionDual for ChainFun fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); - WeightFunctionInfo::new(p.component_index(), true) + WeightFunctionInfo::new(p.component_index.clone(), true) .add( WeightFunction { prefactor: p.m.mapv(|m| m.into()) / (&d * 8.0), diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 79c54864a..e604c3e5f 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -2,10 +2,8 @@ use super::eos::GcPcSaftOptions; use feos_core::parameter::ParameterHetero; use feos_core::MolarWeight; use feos_dft::adsorption::FluidParameters; -use feos_dft::fundamental_measure_theory::{ - FMTContribution, FMTProperties, FMTVersion, MonomerShape, -}; use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; +use feos_saft::{FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; use ndarray::Array1; use num_dual::DualNum; use petgraph::graph::UnGraph; @@ -124,14 +122,10 @@ impl MolarWeight for GcPcSaftFunctional { } } -impl FMTProperties for GcPcSaftFunctionalParameters { - fn component_index(&self) -> Array1 { - self.component_index.clone() - } - +impl HardSphereProperties for GcPcSaftFunctionalParameters { fn monomer_shape>(&self, _: N) -> MonomerShape { let m = self.m.mapv(N::from); - MonomerShape::Heterosegmented([m.clone(), m.clone(), m.clone(), m]) + MonomerShape::Heterosegmented([m.clone(), m.clone(), m.clone(), m], &self.component_index) } fn hs_diameter>(&self, temperature: D) -> Array1 { diff --git a/src/gc_pcsaft/eos/association.rs b/src/gc_pcsaft/eos/association.rs index d20db7935..9f9ef169a 100644 --- a/src/gc_pcsaft/eos/association.rs +++ b/src/gc_pcsaft/eos/association.rs @@ -1,6 +1,6 @@ -use super::hard_sphere::zeta; use super::GcPcSaftEosParameters; use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::linalg::{norm, LU}; use num_dual::*; @@ -48,8 +48,8 @@ impl> HelmholtzEnergyDual for Association { let diameter = p.hs_diameter(state.temperature); // Packing fractions - let n2 = zeta(p, &diameter, &state.partial_density, [2])[0] * 6.0; - let n3 = zeta(p, &diameter, &state.partial_density, [3])[0]; + let [zeta2, n3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); + let n2 = zeta2 * 6.0; let n3i = (-n3 + 1.0).recip(); // association strength @@ -116,8 +116,8 @@ impl + ScalarOperand> HelmholtzEnergyDual for CrossAssociatio let diameter = p.hs_diameter(state.temperature); // Packing fractions - let n2 = zeta(p, &diameter, &state.partial_density, [2])[0] * 6.0; - let n3 = zeta(p, &diameter, &state.partial_density, [3])[0]; + let [zeta2, n3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); + let n2 = zeta2 * 6.0; let n3i = (-n3 + 1.0).recip(); // extract densities of associating segments diff --git a/src/gc_pcsaft/eos/dispersion.rs b/src/gc_pcsaft/eos/dispersion.rs index e0e8280dd..b6db3f5ee 100644 --- a/src/gc_pcsaft/eos/dispersion.rs +++ b/src/gc_pcsaft/eos/dispersion.rs @@ -1,6 +1,6 @@ -use super::hard_sphere::zeta; use super::GcPcSaftEosParameters; use feos_core::{HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use num_dual::DualNum; use std::f64::consts::PI; use std::fmt; @@ -74,8 +74,7 @@ impl> HelmholtzEnergyDual for Dispersion { let rho = &state.partial_density; // packing fraction - let diameter = p.hs_diameter(state.temperature); - let eta = zeta(p, &diameter, &state.partial_density, [3])[0]; + let eta = p.zeta(state.temperature, &state.partial_density, [3])[0]; // mean segment number let m = diff --git a/src/gc_pcsaft/eos/hard_chain.rs b/src/gc_pcsaft/eos/hard_chain.rs index 273a7a1d7..ffde39044 100644 --- a/src/gc_pcsaft/eos/hard_chain.rs +++ b/src/gc_pcsaft/eos/hard_chain.rs @@ -1,6 +1,6 @@ -use super::hard_sphere::zeta; use super::GcPcSaftEosParameters; use feos_core::{HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use num_dual::*; use std::fmt; use std::rc::Rc; @@ -16,7 +16,9 @@ impl> HelmholtzEnergyDual for HardChain { let diameter = self.parameters.hs_diameter(state.temperature); // Packing fractions - let [zeta2, zeta3] = zeta(&self.parameters, &diameter, &state.partial_density, [2, 3]); + let [zeta2, zeta3] = + self.parameters + .zeta(state.temperature, &state.partial_density, [2, 3]); // Helmholtz energy let frac_1mz3 = -(zeta3 - 1.0).recip(); diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index 574bdfa82..ddc21cb6e 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -1,6 +1,7 @@ use feos_core::joback::Joback; use feos_core::parameter::ParameterHetero; use feos_core::{EquationOfState, HelmholtzEnergy, IdealGasContribution, MolarWeight}; +use feos_saft::HardSphere; use ndarray::Array1; use quantity::si::*; use std::f64::consts::FRAC_PI_6; @@ -9,13 +10,11 @@ use std::rc::Rc; pub(crate) mod association; pub(crate) mod dispersion; mod hard_chain; -mod hard_sphere; mod parameter; mod polar; use association::{Association, CrossAssociation}; use dispersion::Dispersion; use hard_chain::HardChain; -use hard_sphere::HardSphere; pub use parameter::{GcPcSaftChemicalRecord, GcPcSaftEosParameters}; use polar::Dipole; diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs index 4d7f2039e..d68c12932 100644 --- a/src/gc_pcsaft/eos/parameter.rs +++ b/src/gc_pcsaft/eos/parameter.rs @@ -4,8 +4,10 @@ use feos_core::parameter::{ BinaryRecord, ChemicalRecord, FromSegments, Identifier, ParameterError, ParameterHetero, SegmentCount, SegmentRecord, }; +use feos_saft::{HardSphereProperties, MonomerShape}; use indexmap::IndexMap; use ndarray::{Array1, Array2}; +use num_dual::DualNum; use quantity::si::{JOULE, KB, KELVIN}; use std::borrow::Cow; use std::collections::HashMap; @@ -305,6 +307,20 @@ impl ParameterHetero for GcPcSaftEosParameters { } } +impl HardSphereProperties for GcPcSaftEosParameters { + fn monomer_shape>(&self, _: N) -> MonomerShape { + let m = self.m.mapv(N::from); + MonomerShape::Heterosegmented([m.clone(), m.clone(), m.clone(), m], &self.component_index) + } + + fn hs_diameter>(&self, temperature: D) -> Array1 { + let ti = temperature.recip() * -3.0; + Array1::from_shape_fn(self.sigma.len(), |i| { + -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] + }) + } +} + impl GcPcSaftEosParameters { pub fn to_markdown(&self) -> String { let mut output = String::new(); diff --git a/src/gc_pcsaft/eos/polar.rs b/src/gc_pcsaft/eos/polar.rs index 26cb01913..cce3d9c3d 100644 --- a/src/gc_pcsaft/eos/polar.rs +++ b/src/gc_pcsaft/eos/polar.rs @@ -1,5 +1,6 @@ use super::GcPcSaftEosParameters; use feos_core::{HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use ndarray::prelude::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; diff --git a/src/pcsaft/dft/association.rs b/src/pcsaft/dft/association.rs index 50fc0a4ea..a5fadf396 100644 --- a/src/pcsaft/dft/association.rs +++ b/src/pcsaft/dft/association.rs @@ -6,6 +6,7 @@ use feos_core::EosError; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::DualNum; use std::f64::consts::PI; diff --git a/src/pcsaft/dft/dispersion.rs b/src/pcsaft/dft/dispersion.rs index afe7b3d7e..eb741da5c 100644 --- a/src/pcsaft/dft/dispersion.rs +++ b/src/pcsaft/dft/dispersion.rs @@ -5,6 +5,7 @@ use feos_core::EosError; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; diff --git a/src/pcsaft/dft/hard_chain.rs b/src/pcsaft/dft/hard_chain.rs index 142447910..f2056b793 100644 --- a/src/pcsaft/dft/hard_chain.rs +++ b/src/pcsaft/dft/hard_chain.rs @@ -1,9 +1,9 @@ use super::PcSaftParameters; use feos_core::EosError; -use feos_dft::fundamental_measure_theory::FMTProperties; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::DualNum; use std::fmt; @@ -24,7 +24,7 @@ impl + ScalarOperand> FunctionalContributionDual for ChainFun fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); - WeightFunctionInfo::new(p.component_index(), true) + WeightFunctionInfo::new(p.component_index().into_owned(), true) .add( WeightFunction { prefactor: p.m.mapv(|m| m.into()) / (&d * 8.0), diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index 43d185aa6..6e949a948 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -6,14 +6,11 @@ use feos_core::joback::Joback; use feos_core::parameter::Parameter; use feos_core::{IdealGasContribution, MolarWeight}; use feos_dft::adsorption::FluidParameters; -use feos_dft::fundamental_measure_theory::{ - FMTContribution, FMTProperties, FMTVersion, MonomerShape, -}; use feos_dft::solvation::PairPotential; use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; +use feos_saft::{FMTContribution, FMTVersion}; use hard_chain::ChainFunctional; -use ndarray::{Array, Array1, Array2}; -use num_dual::DualNum; +use ndarray::{Array1, Array2}; use num_traits::One; use pure_saft_functional::*; use quantity::si::*; @@ -140,20 +137,6 @@ impl MolarWeight for PcSaftFunctional { } } -impl FMTProperties for PcSaftParameters { - fn component_index(&self) -> Array1 { - Array::from_shape_fn(self.m.len(), |i| i) - } - - fn monomer_shape>(&self, _: N) -> MonomerShape { - MonomerShape::NonSpherical(self.m.mapv(N::from)) - } - - fn hs_diameter>(&self, temperature: D) -> Array1 { - self.hs_diameter(temperature) - } -} - impl FluidParameters for PcSaftFunctional { fn epsilon_k_ff(&self) -> Array1 { self.parameters.epsilon_k.clone() diff --git a/src/pcsaft/dft/polar.rs b/src/pcsaft/dft/polar.rs index 2b724f3ac..9534241f6 100644 --- a/src/pcsaft/dft/polar.rs +++ b/src/pcsaft/dft/polar.rs @@ -3,6 +3,7 @@ use crate::pcsaft::eos::polar::{ MeanSegmentNumbers, Multipole, AD, ADQ, ALPHA, AQ, BD, BDQ, BQ, CD, CDQ, CQ, PI_SQ_43, }; use feos_core::EosError; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index 82cc1d65f..d60dd88d0 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -5,10 +5,10 @@ use crate::pcsaft::eos::association::{assoc_site_frac_a, assoc_site_frac_ab}; use crate::pcsaft::eos::dispersion::{A0, A1, A2, B0, B1, B2}; use crate::pcsaft::eos::polar::{AD, AQ, BD, BQ, CD, CQ, PI_SQ_43}; use feos_core::{EosError, EosResult}; -use feos_dft::fundamental_measure_theory::FMTVersion; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::{FMTVersion, HardSphereProperties}; use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; diff --git a/src/pcsaft/eos/association.rs b/src/pcsaft/eos/association.rs index 3e4b0a5de..655c1c369 100644 --- a/src/pcsaft/eos/association.rs +++ b/src/pcsaft/eos/association.rs @@ -1,5 +1,6 @@ use super::PcSaftParameters; use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::linalg::{norm, LU}; use num_dual::*; diff --git a/src/pcsaft/eos/dispersion.rs b/src/pcsaft/eos/dispersion.rs index 24e55e41a..2fe7b8d86 100644 --- a/src/pcsaft/eos/dispersion.rs +++ b/src/pcsaft/eos/dispersion.rs @@ -1,5 +1,6 @@ use super::PcSaftParameters; use feos_core::{HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; use std::fmt; diff --git a/src/pcsaft/eos/hard_chain.rs b/src/pcsaft/eos/hard_chain.rs index 55897c976..f3f5620da 100644 --- a/src/pcsaft/eos/hard_chain.rs +++ b/src/pcsaft/eos/hard_chain.rs @@ -1,6 +1,6 @@ -use super::hard_sphere::zeta; use super::PcSaftParameters; use feos_core::{HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use ndarray::Array; use num_dual::*; use std::fmt; @@ -12,12 +12,13 @@ pub struct HardChain { impl> HelmholtzEnergyDual for HardChain { fn helmholtz_energy(&self, state: &StateHD) -> D { + let p = &self.parameters; let d = self.parameters.hs_diameter(state.temperature); - let zeta = zeta(&self.parameters.m, &state.partial_density, &d); - let frac_1mz3 = -(zeta[3] - 1.0).recip(); - let c = zeta[2] * frac_1mz3 * frac_1mz3; + let [zeta2, zeta3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); + let frac_1mz3 = -(zeta3 - 1.0).recip(); + let c = zeta2 * frac_1mz3 * frac_1mz3; let g_hs = - d.mapv(|d| frac_1mz3 + d * c * 1.5 - d.powi(2) * c.powi(2) * (zeta[3] - 1.0) * 0.5); + d.mapv(|d| frac_1mz3 + d * c * 1.5 - d.powi(2) * c.powi(2) * (zeta3 - 1.0) * 0.5); Array::from_shape_fn(self.parameters.m.len(), |i| { state.partial_density[i] * (1.0 - self.parameters.m[i]) * g_hs[i].ln() }) diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index 13f2d79fd..2af502e7e 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -5,6 +5,7 @@ use feos_core::{ Contributions, EntropyScaling, EosError, EosResult, EquationOfState, HelmholtzEnergy, IdealGasContribution, MolarWeight, State, }; +use feos_saft::HardSphere; use ndarray::Array1; use quantity::si::*; use std::f64::consts::{FRAC_PI_6, PI}; @@ -13,13 +14,11 @@ use std::rc::Rc; pub(crate) mod association; pub(crate) mod dispersion; pub(crate) mod hard_chain; -pub(crate) mod hard_sphere; pub(crate) mod polar; mod qspr; use association::{Association, CrossAssociation}; use dispersion::Dispersion; use hard_chain::HardChain; -use hard_sphere::HardSphere; use polar::{DQVariants, Dipole, DipoleQuadrupole, Quadrupole}; use qspr::QSPR; diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs index fb03c935c..e662ce29e 100644 --- a/src/pcsaft/eos/polar.rs +++ b/src/pcsaft/eos/polar.rs @@ -1,5 +1,6 @@ use super::PcSaftParameters; use feos_core::{HelmholtzEnergyDual, StateHD}; +use feos_saft::HardSphereProperties; use ndarray::prelude::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 9d74c4ff8..3ee93a3c8 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -3,7 +3,9 @@ use feos_core::joback::JobackRecord; use feos_core::parameter::{ FromSegments, FromSegmentsBinary, Parameter, ParameterError, PureRecord, }; +use feos_saft::{HardSphereProperties, MonomerShape}; use ndarray::{Array, Array1, Array2}; +use num_dual::DualNum; use num_traits::Zero; use quantity::si::{JOULE, KB, KELVIN}; use serde::{Deserialize, Serialize}; @@ -507,6 +509,19 @@ impl Parameter for PcSaftParameters { } } +impl HardSphereProperties for PcSaftParameters { + fn monomer_shape>(&self, _: N) -> MonomerShape { + MonomerShape::NonSpherical(self.m.mapv(N::from)) + } + + fn hs_diameter>(&self, temperature: D) -> Array1 { + let ti = temperature.recip() * -3.0; + Array::from_shape_fn(self.sigma.len(), |i| { + -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] + }) + } +} + impl PcSaftParameters { pub fn to_markdown(&self) -> String { let mut output = String::new(); diff --git a/src/python/dft.rs b/src/python/dft.rs index 0bff3d2a7..f6665d3c2 100644 --- a/src/python/dft.rs +++ b/src/python/dft.rs @@ -16,11 +16,11 @@ use crate::pets::python::PyPetsParameters; use crate::pets::{PetsFunctional, PetsOptions}; use feos_core::*; use feos_dft::adsorption::*; -use feos_dft::fundamental_measure_theory::{FMTFunctional, FMTVersion}; use feos_dft::interface::*; use feos_dft::python::*; use feos_dft::solvation::*; use feos_dft::*; +use feos_saft::{FMTFunctional, FMTVersion}; use ndarray::{Array1, Array2}; use numpy::convert::ToPyArray; use numpy::{PyArray1, PyArray2, PyArray4}; diff --git a/tests/pcsaft/dft.rs b/tests/pcsaft/dft.rs index 816c7bb8a..7f21f7a94 100644 --- a/tests/pcsaft/dft.rs +++ b/tests/pcsaft/dft.rs @@ -4,8 +4,8 @@ use approx::assert_relative_eq; use feos::pcsaft::{PcSaft, PcSaftFunctional, PcSaftParameters}; use feos_core::parameter::{IdentifierOption, Parameter}; use feos_core::{Contributions, PhaseEquilibrium, State}; -use feos_dft::fundamental_measure_theory::FMTVersion; use feos_dft::interface::PlanarInterface; +use feos_saft::FMTVersion; use ndarray::{arr1, Axis}; use quantity::si::*; use std::error::Error; From 2631c7c8d5669b51d0d9edfd143e39b31107eb54 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sun, 3 Jul 2022 11:41:10 +0200 Subject: [PATCH 2/9] Use a single generic implementation for association across models --- Cargo.toml | 6 +- feos-saft/Cargo.toml | 4 +- feos-saft/src/fundamental_measure_theory.rs | 346 -------------------- feos-saft/src/hard_sphere.rs | 153 --------- feos-saft/src/lib.rs | 10 +- parameters/pcsaft/sauer2014_hetero.json | 16 +- parameters/pcsaft/sauer2014_homo.json | 36 +- src/gc_pcsaft/dft/association.rs | 198 ----------- src/gc_pcsaft/dft/dispersion.rs | 1 + src/gc_pcsaft/dft/hard_chain.rs | 1 + src/gc_pcsaft/dft/mod.rs | 16 +- src/gc_pcsaft/dft/parameter.rs | 161 ++++----- src/gc_pcsaft/eos/hard_sphere.rs | 65 ++-- src/gc_pcsaft/eos/mod.rs | 27 +- src/gc_pcsaft/eos/parameter.rs | 122 +++---- src/gc_pcsaft/mod.rs | 2 +- src/gc_pcsaft/python/mod.rs | 47 ++- src/gc_pcsaft/record.rs | 45 +-- src/pcsaft/dft/association.rs | 188 ----------- src/pcsaft/dft/mod.rs | 17 +- src/pcsaft/dft/pure_saft_functional.rs | 28 +- src/pcsaft/eos/mod.rs | 24 +- src/pcsaft/eos/polar.rs | 2 + src/pcsaft/eos/qspr.rs | 13 +- src/pcsaft/mod.rs | 4 +- src/pcsaft/parameters.rs | 179 ++++------ src/pcsaft/python.rs | 43 +-- src/python/dft.rs | 12 +- src/python/eos.rs | 12 +- tests/gc_pcsaft/dft.rs | 122 ++++++- tests/pcsaft/test_parameters.json | 8 +- 31 files changed, 495 insertions(+), 1413 deletions(-) delete mode 100644 feos-saft/src/fundamental_measure_theory.rs delete mode 100644 feos-saft/src/hard_sphere.rs delete mode 100644 src/gc_pcsaft/dft/association.rs delete mode 100644 src/pcsaft/dft/association.rs diff --git a/Cargo.toml b/Cargo.toml index e7829b771..f0454b888 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,11 +47,11 @@ optional = true approx = "0.4" [features] -default = [] +default = ["pcsaft", "dft", "gc_pcsaft"] dft = ["feos-dft", "petgraph", "feos-saft?/dft"] estimator = [] -pcsaft = ["feos-saft"] -gc_pcsaft = ["feos-saft"] +pcsaft = ["feos-saft/association"] +gc_pcsaft = ["feos-saft/association"] uvtheory = ["lazy_static"] pets = [] python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python", "feos-saft?/python"] diff --git a/feos-saft/Cargo.toml b/feos-saft/Cargo.toml index 549e80c45..fdc271d69 100644 --- a/feos-saft/Cargo.toml +++ b/feos-saft/Cargo.toml @@ -22,8 +22,10 @@ feos-core = { version = "0.3", path = "../feos-core" } feos-dft = { version = "0.2", path = "../feos-dft", optional = true } ndarray = { version = "0.15", features = ["serde"] } pyo3 = { version = "0.16", optional = true } +serde = "1.0" [features] -default = [] +default = ["association"] +association = [] dft = ["feos-dft"] python = ["pyo3"] \ No newline at end of file diff --git a/feos-saft/src/fundamental_measure_theory.rs b/feos-saft/src/fundamental_measure_theory.rs deleted file mode 100644 index 6e346b810..000000000 --- a/feos-saft/src/fundamental_measure_theory.rs +++ /dev/null @@ -1,346 +0,0 @@ -//! Helmholtz energy functionals from fundamental measure theory. -use feos_core::EosResult; -use feos_dft::adsorption::FluidParameters; -use feos_dft::solvation::PairPotential; -use feos_dft::{ - FunctionalContribution, FunctionalContributionDual, HelmholtzEnergyFunctional, MoleculeShape, - WeightFunction, WeightFunctionInfo, WeightFunctionShape, DFT, -}; -use ndarray::*; -use num_dual::DualNum; -use std::f64::consts::PI; -use std::fmt; -use std::rc::Rc; - -use crate::{HardSphereProperties, MonomerShape}; - -const PI36M1: f64 = 1.0 / (36.0 * PI); -const N3_CUTOFF: f64 = 1e-5; - -/// Different versions of fundamental measure theory -#[derive(Clone, Copy)] -#[cfg_attr(feature = "python", pyo3::pyclass)] -pub enum FMTVersion { - /// White Bear ([Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313)) or modified ([Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)) fundamental measure theory - WhiteBear, - /// Scalar fundamental measure theory by [Kierlik and Rosinberg, 1990](https://doi.org/10.1103/PhysRevA.42.3382) - KierlikRosinberg, - /// Anti-symmetric White Bear fundamental measure theory ([Rosenfeld et al., 1997](https://doi.org/10.1103/PhysRevE.55.4245)) and SI of ([Kessler et al., 2021](https://doi.org/10.1016/j.micromeso.2021.111263)) - AntiSymWhiteBear, -} - -/// The [FunctionalContribution] for the hard sphere functional. -pub struct FMTContribution

{ - pub properties: Rc

, - version: FMTVersion, -} - -impl

Clone for FMTContribution

{ - fn clone(&self) -> Self { - Self { - properties: self.properties.clone(), - version: self.version, - } - } -} - -impl

FMTContribution

{ - pub fn new(properties: &Rc

, version: FMTVersion) -> Self { - Self { - properties: properties.clone(), - version, - } - } -} - -impl> FunctionalContributionDual - for FMTContribution

-{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { - let r = self.properties.hs_diameter(temperature) * 0.5; - let [c0, c1, c2, c3] = self.properties.geometry_coefficients(temperature); - match (self.version, r.len()) { - (FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear, 1) => { - WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) - .extend( - vec![ - WeightFunctionShape::Delta, - WeightFunctionShape::Theta, - WeightFunctionShape::DeltaVec, - ] - .into_iter() - .zip([c2, c3.clone(), c3]) - .map(|(s, c)| WeightFunction { - prefactor: c, - kernel_radius: r.clone(), - shape: s, - }) - .collect(), - false, - ) - } - (FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear, _) => { - WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) - .add( - WeightFunction { - prefactor: Zip::from(&c0) - .and(&r) - .map_collect(|&c, &r| r.powi(-2) * c / (4.0 * PI)), - kernel_radius: r.clone(), - shape: WeightFunctionShape::Delta, - }, - true, - ) - .add( - WeightFunction { - prefactor: Zip::from(&c1) - .and(&r) - .map_collect(|&c, &r| r.recip() * c / (4.0 * PI)), - kernel_radius: r.clone(), - shape: WeightFunctionShape::Delta, - }, - true, - ) - .add( - WeightFunction { - prefactor: c2, - kernel_radius: r.clone(), - shape: WeightFunctionShape::Delta, - }, - true, - ) - .add( - WeightFunction { - prefactor: c3.clone(), - kernel_radius: r.clone(), - shape: WeightFunctionShape::Theta, - }, - true, - ) - .add( - WeightFunction { - prefactor: Zip::from(&c3) - .and(&r) - .map_collect(|&c, &r| r.recip() * c / (4.0 * PI)), - kernel_radius: r.clone(), - shape: WeightFunctionShape::DeltaVec, - }, - true, - ) - .add( - WeightFunction { - prefactor: c3, - kernel_radius: r, - shape: WeightFunctionShape::DeltaVec, - }, - true, - ) - } - (FMTVersion::KierlikRosinberg, _) => { - WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) - .extend( - vec![ - WeightFunctionShape::KR0, - WeightFunctionShape::KR1, - WeightFunctionShape::Delta, - WeightFunctionShape::Theta, - ] - .into_iter() - .zip(self.properties.geometry_coefficients(temperature)) - .map(|(s, c)| WeightFunction { - prefactor: c, - kernel_radius: r.clone(), - shape: s, - }) - .collect(), - true, - ) - } - } - } - - fn calculate_helmholtz_energy_density( - &self, - temperature: N, - weighted_densities: ArrayView2, - ) -> EosResult> { - let pure_component_weighted_densities = matches!( - self.version, - FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear - ) && self.properties.component_index().len() == 1; - - // scalar weighted densities - let (n2, n3) = if pure_component_weighted_densities { - ( - weighted_densities.index_axis(Axis(0), 0), - weighted_densities.index_axis(Axis(0), 1), - ) - } else { - ( - weighted_densities.index_axis(Axis(0), 2), - weighted_densities.index_axis(Axis(0), 3), - ) - }; - - let (n0, n1) = if pure_component_weighted_densities { - let r = self.properties.hs_diameter(temperature)[0] * 0.5; - ( - n2.mapv(|n2| n2 / (r * r * 4.0 * PI)), - n2.mapv(|n2| n2 / (r * 4.0 * PI)), - ) - } else { - ( - weighted_densities.index_axis(Axis(0), 0).to_owned(), - weighted_densities.index_axis(Axis(0), 1).to_owned(), - ) - }; - - // vector weighted densities - let (n1n2, n2n2) = match self.version { - FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear => { - let (n1v, n2v) = if pure_component_weighted_densities { - let r = self.properties.hs_diameter(temperature)[0] * 0.5; - let n2v = weighted_densities.slice_axis(Axis(0), Slice::new(2, None, 1)); - (n2v.mapv(|n2v| n2v / (r * 4.0 * PI)), n2v) - } else { - let dim = ((weighted_densities.shape()[0] - 4) / 2) as isize; - ( - weighted_densities - .slice_axis(Axis(0), Slice::new(4, Some(4 + dim), 1)) - .to_owned(), - weighted_densities - .slice_axis(Axis(0), Slice::new(4 + dim, Some(4 + 2 * dim), 1)), - ) - }; - match self.version { - FMTVersion::WhiteBear => ( - &n1 * &n2 - (&n1v * &n2v).sum_axis(Axis(0)), - &n2 * &n2 - (&n2v * &n2v).sum_axis(Axis(0)) * 3.0, - ), - FMTVersion::AntiSymWhiteBear => { - let mut xi2 = (&n2v * &n2v).sum_axis(Axis(0)) / n2.map(|n| n.powi(2)); - xi2.iter_mut().for_each(|x| { - if x.re() > 1.0 { - *x = N::one() - } - }); - ( - &n1 * &n2 - (&n1v * &n2v).sum_axis(Axis(0)), - &n2 * &n2 * xi2.mapv(|x| (-x + 1.0).powi(3)), - ) - } - _ => unreachable!(), - } - } - FMTVersion::KierlikRosinberg => (&n1 * &n2, &n2 * &n2), - }; - - // auxiliary variables - let ln31 = n3.mapv(|n3| (-n3).ln_1p()); - let n3rec = n3.mapv(|n3| n3.recip()); - let n3m1 = n3.mapv(|n3| -n3 + 1.0); - let n3m1rec = n3m1.mapv(|n3m1| n3m1.recip()); - - // use Taylor expansion for f3 at low densities to avoid numerical issues - let mut f3 = (&n3m1 * &n3m1 * &ln31 + n3) * &n3rec * n3rec * &n3m1rec * &n3m1rec; - f3.iter_mut().zip(n3).for_each(|(f3, &n3)| { - if n3.re() < N3_CUTOFF { - *f3 = (((n3 * 35.0 / 6.0 + 4.8) * n3 + 3.75) * n3 + 8.0 / 3.0) * n3 + 1.5; - } - }); - Ok(-(&n0 * &ln31) + n1n2 * &n3m1rec + n2n2 * n2 * PI36M1 * f3) - } -} - -impl fmt::Display for FMTContribution

{ - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let ver = match self.version { - FMTVersion::WhiteBear => "WB", - FMTVersion::KierlikRosinberg => "KR", - FMTVersion::AntiSymWhiteBear => "AntiSymWB", - }; - write!(f, "FMT functional ({})", ver) - } -} - -struct HardSphereParameters { - sigma: Array1, -} - -impl HardSphereProperties for HardSphereParameters { - fn monomer_shape(&self, _: N) -> MonomerShape { - MonomerShape::Spherical(self.sigma.len()) - } - - fn hs_diameter>(&self, _: N) -> Array1 { - self.sigma.mapv(N::from) - } -} - -/// [HelmholtzEnergyFunctional] for hard sphere systems. -pub struct FMTFunctional { - properties: Rc, - contributions: Vec>, - version: FMTVersion, -} - -impl FMTFunctional { - pub fn new(sigma: &Array1, version: FMTVersion) -> DFT { - let properties = Rc::new(HardSphereParameters { - sigma: sigma.clone(), - }); - let contributions: Vec> = - vec![Box::new(FMTContribution::new(&properties, version))]; - (Self { - properties, - contributions, - version, - }) - .into() - } -} - -impl HelmholtzEnergyFunctional for FMTFunctional { - fn contributions(&self) -> &[Box] { - &self.contributions - } - - fn subset(&self, component_list: &[usize]) -> DFT { - let sigma = component_list - .iter() - .map(|&c| self.properties.sigma[c]) - .collect(); - Self::new(&sigma, self.version) - } - - fn compute_max_density(&self, moles: &Array1) -> f64 { - moles.sum() / (moles * &self.properties.sigma).sum() * 1.2 - } - - fn molecule_shape(&self) -> MoleculeShape { - MoleculeShape::Spherical(self.properties.sigma.len()) - } -} - -impl PairPotential for FMTFunctional { - fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Array2 { - let s = &self.properties.sigma; - Array::from_shape_fn((s.len(), r.len()), |(j, k)| { - if r[k] > 0.5 * (s[i] + s[j]) { - 0.0 - } else { - f64::INFINITY - } - }) - } -} - -impl FluidParameters for FMTFunctional { - fn epsilon_k_ff(&self) -> Array1 { - Array::zeros(self.properties.sigma.len()) - } - - fn sigma_ff(&self) -> &Array1 { - &self.properties.sigma - } -} diff --git a/feos-saft/src/hard_sphere.rs b/feos-saft/src/hard_sphere.rs deleted file mode 100644 index 5bf17e16c..000000000 --- a/feos-saft/src/hard_sphere.rs +++ /dev/null @@ -1,153 +0,0 @@ -use feos_core::{HelmholtzEnergyDual, StateHD}; -use ndarray::*; -use num_dual::DualNum; -use std::borrow::Cow; -use std::f64::consts::FRAC_PI_6; -use std::fmt; -use std::rc::Rc; - -/// Different monomer shapes for FMT. -pub enum MonomerShape<'a, D> { - /// For spherical monomers, the number of components. - Spherical(usize), - /// For non-spherical molecules in a homosegmented approach, the - /// chain length parameter $m$. - NonSpherical(Array1), - /// For non-spherical molecules in a heterosegmented approach, - /// the geometry factors for every segment and the component - /// index for every segment. - Heterosegmented([Array1; 4], &'a Array1), -} - -/// Properties of (generalized) hard sphere systems. -pub trait HardSphereProperties { - fn monomer_shape>(&self, temperature: D) -> MonomerShape; - fn hs_diameter>(&self, temperature: D) -> Array1; - - fn component_index(&self) -> Cow> { - match self.monomer_shape(1.0) { - MonomerShape::Spherical(n) => Cow::Owned(Array1::from_shape_fn(n, |i| i)), - MonomerShape::NonSpherical(m) => Cow::Owned(Array1::from_shape_fn(m.len(), |i| i)), - MonomerShape::Heterosegmented(_, component_index) => Cow::Borrowed(component_index), - } - } - - fn geometry_coefficients>(&self, temperature: D) -> [Array1; 4] { - match self.monomer_shape(temperature) { - MonomerShape::Spherical(n) => { - let m = Array1::ones(n); - [m.clone(), m.clone(), m.clone(), m] - } - MonomerShape::NonSpherical(m) => [m.clone(), m.clone(), m.clone(), m], - MonomerShape::Heterosegmented(g, _) => g, - } - } - - fn zeta, const N: usize>( - &self, - temperature: D, - partial_density: &Array1, - k: [i32; N], - ) -> [D; N] { - let component_index = self.component_index(); - let geometry_coefficients = self.geometry_coefficients(temperature); - let diameter = self.hs_diameter(temperature); - let mut zeta = [D::zero(); N]; - for i in 0..diameter.len() { - for (z, &k) in zeta.iter_mut().zip(k.iter()) { - *z += partial_density[component_index[i]] - * diameter[i].powi(k) - * (geometry_coefficients[k as usize][i] * FRAC_PI_6); - } - } - - zeta - } - - fn zeta_23>(&self, temperature: D, molefracs: &Array1) -> D { - let component_index = self.component_index(); - let geometry_coefficients = self.geometry_coefficients(temperature); - let diameter = self.hs_diameter(temperature); - let mut zeta: [D; 2] = [D::zero(); 2]; - for i in 0..diameter.len() { - for (k, z) in zeta.iter_mut().enumerate() { - *z += molefracs[component_index[i]] - * diameter[i].powi((k + 2) as i32) - * (geometry_coefficients[k + 2][i] * FRAC_PI_6); - } - } - - zeta[0] / zeta[1] - } -} - -pub struct HardSphere

{ - pub parameters: Rc

, -} - -impl, P: HardSphereProperties> HelmholtzEnergyDual for HardSphere

{ - fn helmholtz_energy(&self, state: &StateHD) -> D { - let p = &self.parameters; - let zeta = p.zeta(state.temperature, &state.partial_density, [0, 1, 2, 3]); - let frac_1mz3 = -(zeta[3] - 1.0).recip(); - let zeta_23 = p.zeta_23(state.temperature, &state.molefracs); - state.volume * 6.0 / std::f64::consts::PI - * (zeta[1] * zeta[2] * frac_1mz3 * 3.0 - + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23 - + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p()) - } -} - -impl

fmt::Display for HardSphere

{ - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Hard Sphere") - } -} - -// #[cfg(test)] -// mod tests { -// use super::*; -// use crate::pcsaft::parameters::utils::{ -// butane_parameters, propane_butane_parameters, propane_parameters, -// }; -// use approx::assert_relative_eq; -// use ndarray::arr1; - -// #[test] -// fn helmholtz_energy() { -// let hs = HardSphere { -// parameters: propane_parameters(), -// }; -// let t = 250.0; -// let v = 1000.0; -// let n = 1.0; -// let s = StateHD::new(t, v, arr1(&[n])); -// let a_rust = hs.helmholtz_energy(&s); -// assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10); -// } - -// #[test] -// fn mix() { -// let c1 = HardSphere { -// parameters: propane_parameters(), -// }; -// let c2 = HardSphere { -// parameters: butane_parameters(), -// }; -// let c12 = HardSphere { -// parameters: propane_butane_parameters(), -// }; -// let t = 250.0; -// let v = 2.5e28; -// let n = 1.0; -// let s = StateHD::new(t, v, arr1(&[n])); -// let a1 = c1.helmholtz_energy(&s); -// let a2 = c2.helmholtz_energy(&s); -// let s1m = StateHD::new(t, v, arr1(&[n, 0.0])); -// let a1m = c12.helmholtz_energy(&s1m); -// let s2m = StateHD::new(t, v, arr1(&[0.0, n])); -// let a2m = c12.helmholtz_energy(&s2m); -// assert_relative_eq!(a1, a1m, epsilon = 1e-14); -// assert_relative_eq!(a2, a2m, epsilon = 1e-14); -// } -// } diff --git a/feos-saft/src/lib.rs b/feos-saft/src/lib.rs index ce09bc930..fe2d0d59f 100644 --- a/feos-saft/src/lib.rs +++ b/feos-saft/src/lib.rs @@ -1,7 +1,9 @@ mod hard_sphere; +#[cfg(feature = "dft")] +pub use hard_sphere::{FMTContribution, FMTFunctional, FMTVersion}; pub use hard_sphere::{HardSphere, HardSphereProperties, MonomerShape}; -#[cfg(feature = "dft")] -mod fundamental_measure_theory; -#[cfg(feature = "dft")] -pub use fundamental_measure_theory::{FMTContribution, FMTFunctional, FMTVersion}; +#[cfg(feature = "association")] +mod association; +#[cfg(feature = "association")] +pub use association::{Association, AssociationParameters, AssociationRecord}; diff --git a/parameters/pcsaft/sauer2014_hetero.json b/parameters/pcsaft/sauer2014_hetero.json index da9755997..d5c0aa092 100644 --- a/parameters/pcsaft/sauer2014_hetero.json +++ b/parameters/pcsaft/sauer2014_hetero.json @@ -191,8 +191,12 @@ "m": 1.0231, "sigma": 2.7702, "epsilon_k": 334.29, - "epsilon_k_ab": 2575.9, - "kappa_ab": 0.009583 + "association_record": { + "epsilon_k_ab": 2575.9, + "kappa_ab": 0.009583, + "na": 1.0, + "nb": 1.0 + } }, "molarweight": 17.00734 }, @@ -202,8 +206,12 @@ "m": 0.82284, "sigma": 3.1129, "epsilon_k": 309.93, - "epsilon_k_ab": 1471.5, - "kappa_ab": 0.005769 + "association_record": { + "epsilon_k_ab": 1471.5, + "kappa_ab": 0.005769, + "na": 1.0, + "nb": 1.0 + } }, "molarweight": 16.02238 } diff --git a/parameters/pcsaft/sauer2014_homo.json b/parameters/pcsaft/sauer2014_homo.json index 67aa45680..534105508 100644 --- a/parameters/pcsaft/sauer2014_homo.json +++ b/parameters/pcsaft/sauer2014_homo.json @@ -17,7 +17,6 @@ }, "molarweight": 14.02658 }, - { "identifier": ">CH", "model_record": { @@ -57,7 +56,7 @@ { "identifier": "=C<", "model_record": { - "m": 0.86367, + "m": 0.86367, "sigma": 3.1815, "epsilon_k": 156.31 }, @@ -66,7 +65,7 @@ { "identifier": "C≡CH", "model_record": { - "m": 1.3279, + "m": 1.3279, "sigma": 2.9421, "epsilon_k": 223.05 }, @@ -75,7 +74,7 @@ { "identifier": "CH2_hex", "model_record": { - "m": 0.39496, + "m": 0.39496, "sigma": 3.9126, "epsilon_k": 289.03 }, @@ -84,7 +83,7 @@ { "identifier": "CH_hex", "model_record": { - "m": 0.02880, + "m": 0.02880, "sigma": 8.9779, "epsilon_k": 1306.7 }, @@ -93,7 +92,7 @@ { "identifier": "CH2_pent", "model_record": { - "m": 0.46742, + "m": 0.46742, "sigma": 3.7272, "epsilon_k": 267.16 }, @@ -102,7 +101,7 @@ { "identifier": "CH_pent", "model_record": { - "m": 0.03314, + "m": 0.03314, "sigma": 7.7190, "epsilon_k": 1297.7 }, @@ -111,7 +110,7 @@ { "identifier": "CH_arom", "model_record": { - "m": 0.42335, + "m": 0.42335, "sigma": 3.7270, "epsilon_k": 274.41 }, @@ -169,7 +168,7 @@ { "identifier": "HCOO", "model_record": { - "m": 1.7525, + "m": 1.7525, "sigma": 2.9043, "epsilon_k": 229.63, "mu": 2.7916 @@ -192,8 +191,12 @@ "m": 0.40200, "sigma": 3.2859, "epsilon_k": 488.66, - "epsilon_k_ab": 2517.0, - "kappa_ab": 0.006825 + "association_record": { + "epsilon_k_ab": 2517.0, + "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0 + } }, "molarweight": 17.00734 }, @@ -203,10 +206,13 @@ "m": 0.40558, "sigma": 3.6456, "epsilon_k": 467.59, - "epsilon_k_ab": 1064.6, - "kappa_ab": 0.026662 + "association_record": { + "epsilon_k_ab": 1064.6, + "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0 + } }, "molarweight": 16.02238 } -] - +] \ No newline at end of file diff --git a/src/gc_pcsaft/dft/association.rs b/src/gc_pcsaft/dft/association.rs deleted file mode 100644 index 540414e58..000000000 --- a/src/gc_pcsaft/dft/association.rs +++ /dev/null @@ -1,198 +0,0 @@ -use super::parameter::GcPcSaftFunctionalParameters; -use crate::gc_pcsaft::eos::association::{ - assoc_site_frac_a, assoc_site_frac_ab, helmholtz_energy_density_cross_association, -}; -use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; -use ndarray::*; -use num_dual::DualNum; -use std::f64::consts::PI; -use std::fmt; -use std::ops::MulAssign; -use std::rc::Rc; - -pub const N0_CUTOFF: f64 = 1e-9; - -#[derive(Clone)] -pub struct AssociationFunctional { - parameters: Rc, - max_iter: usize, - tol: f64, -} - -impl AssociationFunctional { - pub fn new(parameters: &Rc, max_iter: usize, tol: f64) -> Self { - Self { - parameters: parameters.clone(), - max_iter, - tol, - } - } -} - -impl FunctionalContributionDual for AssociationFunctional -where - N: DualNum + ScalarOperand, -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { - let p = &self.parameters; - let r = p.hs_diameter(temperature) * 0.5; - WeightFunctionInfo::new(p.component_index.clone(), false) - .add( - WeightFunction::new_scaled(r.clone(), WeightFunctionShape::Delta), - false, - ) - .add( - WeightFunction { - prefactor: p.m.mapv(N::from), - kernel_radius: r.clone(), - shape: WeightFunctionShape::DeltaVec, - }, - false, - ) - .add( - WeightFunction { - prefactor: p.m.mapv(N::from), - kernel_radius: r, - shape: WeightFunctionShape::Theta, - }, - true, - ) - } - - fn calculate_helmholtz_energy_density( - &self, - temperature: N, - weighted_densities: ArrayView2, - ) -> Result, EosError> { - let p = &self.parameters; - - // number of segments - let segments = p.m.len(); - - // number of associating segments - let nassoc = p.assoc_segment.len(); - - // number of dimensions - let dim = (weighted_densities.shape()[0] - 1) / segments - 1; - - // weighted densities - let n0i = weighted_densities.slice_axis(Axis(0), Slice::new(0, Some(segments as isize), 1)); - let n2vi: Vec<_> = (0..dim) - .map(|i| { - weighted_densities.slice_axis( - Axis(0), - Slice::new( - (segments * (i + 1)) as isize, - Some((segments * (i + 2)) as isize), - 1, - ), - ) - }) - .collect(); - let n3 = weighted_densities.index_axis(Axis(0), segments * (dim + 1)); - - // calculate rho0 (only associating segments) - let diameter = p.hs_diameter(temperature); - let mut n2i = n0i.to_owned(); - for (i, mut n2i) in n2i.outer_iter_mut().enumerate() { - n2i.mul_assign(diameter[i].powi(2) * p.m[i] * PI); - } - let mut rho0: Array2 = (n2vi - .iter() - .fold(Array2::zeros(n2i.raw_dim()), |acc, n2vi| acc + n2vi * n2vi) - / -(&n2i * &n2i) - + 1.0) - * n0i; - rho0.iter_mut().zip(&n0i).for_each(|(rho0, &n0i)| { - if n0i.re() < N0_CUTOFF { - *rho0 = n0i; - } - }); - let rho0 = - Array2::from_shape_fn((nassoc, n3.len()), |(i, j)| rho0[(p.assoc_segment[i], j)]); - - // calculate xi - let n2v: Vec<_> = n2vi.iter().map(|n2vi| n2vi.sum_axis(Axis(0))).collect(); - let n2 = n2i.sum_axis(Axis(0)); - let mut xi = n2v - .iter() - .fold(Array::zeros(n2.raw_dim()), |acc, n2v| acc + n2v * n2v) - / -(&n2 * &n2) - + 1.0; - xi.iter_mut() - .zip(&n0i.sum_axis(Axis(0))) - .for_each(|(xi, &n0i)| { - if n0i.re() < N0_CUTOFF { - *xi = N::one(); - } - }); - - // auxiliary variables - let n3i = n3.mapv(|n3| (-n3 + 1.0).recip()); - - // only one associating component - if nassoc == 1 { - // association strength - let k = &n2 * &n3i * diameter[p.assoc_segment[0]] * 0.5; - let deltarho = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) - * ((temperature.recip() * p.epsilon_k_aibj[(0, 0)]).exp_m1() - * p.sigma3_kappa_aibj[(0, 0)]) - * rho0.index_axis(Axis(0), 0); - - let na = p.na[0]; - let nb = p.nb[0]; - let f = |x: N| x.ln() - x * 0.5 + 0.5; - if nb > 0.0 { - // no cross association, two association sites - let xa = deltarho.mapv(|d| assoc_site_frac_ab(d, na, nb)); - let xb = (&xa - 1.0) * (na / nb) + 1.0; - Ok((xa.mapv(f) * na + xb.mapv(f) * nb) * rho0.index_axis(Axis(0), 0)) - } else { - // no cross association, one association site - let xa = deltarho.mapv(|d| assoc_site_frac_a(d, na)); - - Ok(xa.mapv(f) * na * rho0.index_axis(Axis(0), 0)) - } - } else { - let mut x: Array1 = Array::from_elem(2 * nassoc, 0.2); - Ok(rho0 - .view() - .into_shape([nassoc, rho0.len() / nassoc]) - .unwrap() - .axis_iter(Axis(1)) - .zip(n2.iter()) - .zip(n3i.iter()) - .zip(xi.iter()) - .map(|(((rho0, &n2), &n3i), &xi)| { - helmholtz_energy_density_cross_association( - &p.assoc_segment, - &p.sigma3_kappa_aibj, - &p.epsilon_k_aibj, - &p.na, - &p.nb, - temperature, - &rho0, - &diameter, - n2, - n3i, - xi, - self.max_iter, - self.tol, - Some(&mut x), - ) - }) - .collect::, _>>()? - .into_shape(n2.raw_dim()) - .unwrap()) - } - } -} - -impl fmt::Display for AssociationFunctional { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Association functional") - } -} diff --git a/src/gc_pcsaft/dft/dispersion.rs b/src/gc_pcsaft/dft/dispersion.rs index f049ae4bc..f61c5372d 100644 --- a/src/gc_pcsaft/dft/dispersion.rs +++ b/src/gc_pcsaft/dft/dispersion.rs @@ -4,6 +4,7 @@ use feos_core::EosError; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_6, PI}; diff --git a/src/gc_pcsaft/dft/hard_chain.rs b/src/gc_pcsaft/dft/hard_chain.rs index 0b1628489..cf0c5ddb9 100644 --- a/src/gc_pcsaft/dft/hard_chain.rs +++ b/src/gc_pcsaft/dft/hard_chain.rs @@ -3,6 +3,7 @@ use feos_core::EosError; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::HardSphereProperties; use ndarray::*; use num_dual::DualNum; use petgraph::visit::EdgeRef; diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index e604c3e5f..65c999097 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -3,7 +3,7 @@ use feos_core::parameter::ParameterHetero; use feos_core::MolarWeight; use feos_dft::adsorption::FluidParameters; use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use feos_saft::{FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; +use feos_saft::{Association, FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; use ndarray::Array1; use num_dual::DualNum; use petgraph::graph::UnGraph; @@ -11,11 +11,11 @@ use quantity::si::{SIArray1, SIUnit, GRAM, MOL}; use std::f64::consts::FRAC_PI_6; use std::rc::Rc; -mod association; +// mod association; mod dispersion; mod hard_chain; mod parameter; -use association::AssociationFunctional; +// use association::AssociationFunctional; use dispersion::AttractiveFunctional; use hard_chain::ChainFunctional; pub use parameter::GcPcSaftFunctionalParameters; @@ -57,9 +57,10 @@ impl GcPcSaftFunctional { contributions.push(Box::new(att)); // Association - if !parameters.assoc_segment.is_empty() { - let assoc = AssociationFunctional::new( + if !parameters.association.assoc_comp.is_empty() { + let assoc = Association::new( ¶meters, + ¶meters.association, saft_options.max_iter_cross_assoc, saft_options.tol_cross_assoc, ); @@ -129,7 +130,10 @@ impl HardSphereProperties for GcPcSaftFunctionalParameters { } fn hs_diameter>(&self, temperature: D) -> Array1 { - self.hs_diameter(temperature) + let ti = temperature.recip() * -3.0; + Array1::from_shape_fn(self.sigma.len(), |i| { + -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] + }) } } diff --git a/src/gc_pcsaft/dft/parameter.rs b/src/gc_pcsaft/dft/parameter.rs index fb65b9bde..573d283c4 100644 --- a/src/gc_pcsaft/dft/parameter.rs +++ b/src/gc_pcsaft/dft/parameter.rs @@ -3,9 +3,9 @@ use feos_core::joback::JobackRecord; use feos_core::parameter::{ BinaryRecord, ChemicalRecord, ParameterError, ParameterHetero, SegmentRecord, }; +use feos_saft::AssociationParameters; use indexmap::IndexMap; use ndarray::{Array1, Array2}; -use num_dual::DualNum; use petgraph::dot::{Config, Dot}; use petgraph::graph::{Graph, UnGraph}; use std::fmt::Write; @@ -22,17 +22,11 @@ pub struct GcPcSaftFunctionalParameters { pub sigma: Array1, pub epsilon_k: Array1, pub bonds: UnGraph<(), ()>, - pub assoc_segment: Array1, - kappa_ab: Array1, - epsilon_k_ab: Array1, - pub na: Array1, - pub nb: Array1, + pub association: AssociationParameters, pub psi_dft: Array1, pub k_ij: Array2, pub sigma_ij: Array2, pub epsilon_k_ij: Array2, - pub sigma3_kappa_aibj: Array2, - pub epsilon_k_aibj: Array2, chemical_records: Vec, segment_records: Vec>, binary_segment_records: Option>>, @@ -63,17 +57,11 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { let mut sigma = Vec::new(); let mut epsilon_k = Vec::new(); let mut bonds = Graph::default(); - let mut assoc_segment = Vec::new(); - let mut kappa_ab = Vec::new(); - let mut epsilon_k_ab = Vec::new(); - let mut na = Vec::new(); - let mut nb = Vec::new(); + let mut association_records = Vec::new(); let mut psi_dft = Vec::new(); let mut segment_index = 0; for (i, chemical_record) in chemical_records.iter().enumerate() { - // let (segment_list, bond_list) = chemical_record.segment_and_bond_list()?; - bonds.extend_with_edges(chemical_record.bonds.iter().map(|x| { ( (segment_index + x[0]) as u32, @@ -93,16 +81,7 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { sigma.push(segment.model_record.sigma); epsilon_k.push(segment.model_record.epsilon_k); - if let (Some(k), Some(e)) = ( - segment.model_record.kappa_ab, - segment.model_record.epsilon_k_ab, - ) { - assoc_segment.push(segment_index); - kappa_ab.push(k); - epsilon_k_ab.push(e); - na.push(segment.model_record.na.unwrap_or(1.0)); - nb.push(segment.model_record.nb.unwrap_or(1.0)); - } + association_records.push(segment.model_record.association_record.clone()); psi_dft.push(segment.model_record.psi_dft.unwrap_or(PSI_GC_DFT)); @@ -143,33 +122,24 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { }); // Association - let sigma3_kappa_aibj = Array2::from_shape_fn([kappa_ab.len(); 2], |(i, j)| { - (sigma[assoc_segment[i]] * sigma[assoc_segment[j]]).powf(1.5) - * (kappa_ab[i] * kappa_ab[j]).sqrt() - }); - let epsilon_k_aibj = Array2::from_shape_fn([epsilon_k_ab.len(); 2], |(i, j)| { - 0.5 * (epsilon_k_ab[i] + epsilon_k_ab[j]) - }); + let sigma = Array1::from_vec(sigma); + let component_index = Array1::from_vec(component_index); + let association = + AssociationParameters::new(&association_records, &sigma, Some(&component_index)); Ok(Self { molarweight, - component_index: Array1::from_vec(component_index), + component_index, identifiers, m: Array1::from_vec(m), - sigma: Array1::from_vec(sigma), + sigma, epsilon_k: Array1::from_vec(epsilon_k), bonds, - assoc_segment: Array1::from_vec(assoc_segment), - kappa_ab: Array1::from_vec(kappa_ab), - epsilon_k_ab: Array1::from_vec(epsilon_k_ab), - na: Array1::from_vec(na), - nb: Array1::from_vec(nb), + association, psi_dft: Array1::from_vec(psi_dft), k_ij, sigma_ij, epsilon_k_ij, - sigma3_kappa_aibj, - epsilon_k_aibj, chemical_records, segment_records, binary_segment_records, @@ -192,57 +162,50 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { } impl GcPcSaftFunctionalParameters { - pub fn hs_diameter>(&self, temperature: D) -> Array1 { - let ti = temperature.recip() * -3.0; - Array1::from_shape_fn(self.sigma.len(), |i| { - -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] - }) - } - - pub fn to_markdown(&self) -> String { - let mut output = String::new(); - let o = &mut output; - write!( - o, - "|component|molarweight|segment|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|" - ) - .unwrap(); - for i in 0..self.m.len() { - let component = if i > 0 && self.component_index[i] == self.component_index[i - 1] { - "|".to_string() - } else { - let pure = &self.chemical_records[self.component_index[i]].identifier; - format!( - "{}|{}", - pure.name - .as_ref() - .unwrap_or(&format!("Component {}", i + 1)), - self.molarweight[self.component_index[i]] - ) - }; - let association = if let Some(a) = self.assoc_segment.iter().position(|&a| a == i) { - format!( - "{}|{}|{}|{}", - self.kappa_ab[a], self.epsilon_k_ab[a], self.na[a], self.nb[a] - ) - } else { - "|||".to_string() - }; - write!( - o, - "\n|{}|{}|{}|{}|{}|{}|||", - component, - self.identifiers[i], - self.m[i], - self.sigma[i], - self.epsilon_k[i], - association - ) - .unwrap(); - } - - output - } + // pub fn to_markdown(&self) -> String { + // let mut output = String::new(); + // let o = &mut output; + // write!( + // o, + // "|component|molarweight|segment|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|" + // ) + // .unwrap(); + // for i in 0..self.m.len() { + // let component = if i > 0 && self.component_index[i] == self.component_index[i - 1] { + // "|".to_string() + // } else { + // let pure = &self.chemical_records[self.component_index[i]].identifier; + // format!( + // "{}|{}", + // pure.name + // .as_ref() + // .unwrap_or(&format!("Component {}", i + 1)), + // self.molarweight[self.component_index[i]] + // ) + // }; + // let association = if let Some(a) = self.assoc_segment.iter().position(|&a| a == i) { + // format!( + // "{}|{}|{}|{}", + // self.kappa_ab[a], self.epsilon_k_ab[a], self.na[a], self.nb[a] + // ) + // } else { + // "|||".to_string() + // }; + // write!( + // o, + // "\n|{}|{}|{}|{}|{}|{}|||", + // component, + // self.identifiers[i], + // self.m[i], + // self.sigma[i], + // self.epsilon_k[i], + // association + // ) + // .unwrap(); + // } + + // output + // } pub fn graph(&self) -> String { let graph = self @@ -261,13 +224,13 @@ impl std::fmt::Display for GcPcSaftFunctionalParameters { write!(f, "\n\tsigma={}", self.sigma)?; write!(f, "\n\tepsilon_k={}", self.epsilon_k)?; write!(f, "\n\tbonds={:?}", self.bonds)?; - if !self.assoc_segment.is_empty() { - write!(f, "\n\tassoc_segment={}", self.assoc_segment)?; - write!(f, "\n\tkappa_ab={}", self.kappa_ab)?; - write!(f, "\n\tepsilon_k_ab={}", self.epsilon_k_ab)?; - write!(f, "\n\tna={}", self.na)?; - write!(f, "\n\tnb={}", self.nb)?; - } + // if !self.assoc_segment.is_empty() { + // write!(f, "\n\tassoc_segment={}", self.assoc_segment)?; + // write!(f, "\n\tkappa_ab={}", self.kappa_ab)?; + // write!(f, "\n\tepsilon_k_ab={}", self.epsilon_k_ab)?; + // write!(f, "\n\tna={}", self.na)?; + // write!(f, "\n\tnb={}", self.nb)?; + // } write!(f, "\n)") } } diff --git a/src/gc_pcsaft/eos/hard_sphere.rs b/src/gc_pcsaft/eos/hard_sphere.rs index 9dd9568cf..c0455c25c 100644 --- a/src/gc_pcsaft/eos/hard_sphere.rs +++ b/src/gc_pcsaft/eos/hard_sphere.rs @@ -23,14 +23,11 @@ pub struct HardSphere { impl> HelmholtzEnergyDual for HardSphere { fn helmholtz_energy(&self, state: &StateHD) -> D { let diameter = self.parameters.hs_diameter(state.temperature); - let zeta = zeta( - &self.parameters, - &diameter, - &state.partial_density, - [0, 1, 2, 3], - ); + let zeta = self + .parameters + .zeta(&diameter, &state.partial_density, [0, 1, 2, 3]); let frac_1mz3 = -(zeta[3] - 1.0).recip(); - let zeta_23 = zeta_23(&self.parameters, &diameter, &state.molefracs); + let zeta_23 = self.parameters.zeta_23(&diameter, &state.molefracs); state.volume * 6.0 / std::f64::consts::PI * (zeta[1] * zeta[2] * frac_1mz3 * 3.0 + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23 @@ -44,39 +41,37 @@ impl fmt::Display for HardSphere { } } -pub(crate) fn zeta, const N: usize>( - parameters: &GcPcSaftEosParameters, - diameter: &Array1, - partial_density: &Array1, - k: [i32; N], -) -> [D; N] { - let mut zeta = [D::zero(); N]; - for i in 0..parameters.m.len() { - for (z, &k) in zeta.iter_mut().zip(k.iter()) { - *z += partial_density[parameters.component_index[i]] - * diameter[i].powi(k) - * (FRAC_PI_6 * parameters.m[i]); +impl GcPcSaftEosParameters { + pub(crate) fn zeta, const N: usize>( + &self, + diameter: &Array1, + partial_density: &Array1, + k: [i32; N], + ) -> [D; N] { + let mut zeta = [D::zero(); N]; + for i in 0..self.m.len() { + for (z, &k) in zeta.iter_mut().zip(k.iter()) { + *z += partial_density[self.component_index[i]] + * diameter[i].powi(k) + * (FRAC_PI_6 * self.m[i]); + } } - } - zeta -} + zeta + } -fn zeta_23>( - parameters: &GcPcSaftEosParameters, - diameter: &Array1, - molefracs: &Array1, -) -> D { - let mut zeta: [D; 2] = [D::zero(); 2]; - for i in 0..parameters.m.len() { - for (k, z) in zeta.iter_mut().enumerate() { - *z += molefracs[parameters.component_index[i]] - * diameter[i].powi((k + 2) as i32) - * (FRAC_PI_6 * parameters.m[i]); + fn zeta_23>(&self, diameter: &Array1, molefracs: &Array1) -> D { + let mut zeta: [D; 2] = [D::zero(); 2]; + for i in 0..self.m.len() { + for (k, z) in zeta.iter_mut().enumerate() { + *z += molefracs[self.component_index[i]] + * diameter[i].powi((k + 2) as i32) + * (FRAC_PI_6 * self.m[i]); + } } - } - zeta[0] / zeta[1] + zeta[0] / zeta[1] + } } #[cfg(test)] diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index ddc21cb6e..09cb0fcd2 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -1,18 +1,18 @@ use feos_core::joback::Joback; use feos_core::parameter::ParameterHetero; use feos_core::{EquationOfState, HelmholtzEnergy, IdealGasContribution, MolarWeight}; -use feos_saft::HardSphere; +use feos_saft::{Association, HardSphere}; use ndarray::Array1; use quantity::si::*; use std::f64::consts::FRAC_PI_6; use std::rc::Rc; -pub(crate) mod association; +// pub(crate) mod association; pub(crate) mod dispersion; mod hard_chain; -mod parameter; +pub(crate) mod parameter; mod polar; -use association::{Association, CrossAssociation}; +// use association::{Association, CrossAssociation}; use dispersion::Dispersion; use hard_chain::HardChain; pub use parameter::{GcPcSaftChemicalRecord, GcPcSaftEosParameters}; @@ -63,17 +63,14 @@ impl GcPcSaft { contributions.push(Box::new(Dispersion { parameters: parameters.clone(), })); - match parameters.assoc_segment.len() { - 0 => (), - 1 => contributions.push(Box::new(Association { - parameters: parameters.clone(), - })), - _ => contributions.push(Box::new(CrossAssociation { - parameters: parameters.clone(), - max_iter: options.max_iter_cross_assoc, - tol: options.tol_cross_assoc, - })), - }; + if !parameters.association.assoc_comp.is_empty() { + contributions.push(Box::new(Association::new( + ¶meters, + ¶meters.association, + options.max_iter_cross_assoc, + options.tol_cross_assoc, + ))); + } if !parameters.dipole_comp.is_empty() { contributions.push(Box::new(Dipole::new(¶meters))) } diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs index d68c12932..ae1d7e66c 100644 --- a/src/gc_pcsaft/eos/parameter.rs +++ b/src/gc_pcsaft/eos/parameter.rs @@ -4,7 +4,7 @@ use feos_core::parameter::{ BinaryRecord, ChemicalRecord, FromSegments, Identifier, ParameterError, ParameterHetero, SegmentCount, SegmentRecord, }; -use feos_saft::{HardSphereProperties, MonomerShape}; +use feos_saft::{AssociationParameters, HardSphereProperties, MonomerShape}; use indexmap::IndexMap; use ndarray::{Array1, Array2}; use num_dual::DualNum; @@ -67,12 +67,7 @@ pub struct GcPcSaftEosParameters { pub epsilon_k: Array1, pub bonds: IndexMap<[usize; 2], f64>, - pub assoc_segment: Array1, - pub n: Array1, - kappa_ab: Array1, - epsilon_k_ab: Array1, - pub na: Array1, - pub nb: Array1, + pub association: AssociationParameters, pub dipole_comp: Array1, mu: Array1, @@ -84,8 +79,6 @@ pub struct GcPcSaftEosParameters { pub k_ij: Array2, pub sigma_ij: Array2, pub epsilon_k_ij: Array2, - pub sigma3_kappa_aibj: Array2, - pub epsilon_k_aibj: Array2, pub chemical_records: Vec, segment_records: Vec>, @@ -113,12 +106,7 @@ impl ParameterHetero for GcPcSaftEosParameters { let mut sigma = Vec::new(); let mut epsilon_k = Vec::new(); let mut bonds = IndexMap::with_capacity(segment_records.len()); - let mut assoc_segment = Vec::new(); - let mut n = Vec::new(); - let mut kappa_ab = Vec::new(); - let mut epsilon_k_ab = Vec::new(); - let mut na = Vec::new(); - let mut nb = Vec::new(); + let mut association_records = Vec::new(); let mut dipole_comp = Vec::new(); let mut mu = Vec::new(); @@ -132,16 +120,6 @@ impl ParameterHetero for GcPcSaftEosParameters { for (i, chemical_record) in chemical_records.iter().cloned().enumerate() { let mut segment_indices = IndexMap::with_capacity(segment_records.len()); let segment_map = chemical_record.segment_map(&segment_records)?; - // let count: IndexMap<_, _> = chemical_record - // .segments - // .iter() - // .map(|(id, &count)| { - // let segment = segment_map - // .get(id) - // .ok_or_else(|| ParameterError::ComponentsNotFound(id.clone()))?; - // Ok((segment, count)) - // }) - // .collect::>()?; let mut m_i = 0.0; let mut sigma_i = 0.0; @@ -159,17 +137,12 @@ impl ParameterHetero for GcPcSaftEosParameters { sigma.push(segment.model_record.sigma); epsilon_k.push(segment.model_record.epsilon_k); - if let (Some(k), Some(e)) = ( - segment.model_record.kappa_ab, - segment.model_record.epsilon_k_ab, - ) { - assoc_segment.push(m.len() - 1); - n.push(count); - kappa_ab.push(k); - epsilon_k_ab.push(e); - na.push(segment.model_record.na.unwrap_or(1.0)); - nb.push(segment.model_record.nb.unwrap_or(1.0)); - } + let mut assoc = segment.model_record.association_record.clone(); + if let Some(mut assoc) = assoc.as_mut() { + assoc.na *= count; + assoc.nb *= count; + }; + association_records.push(assoc); m_i += segment.model_record.m * count; sigma_i += segment.model_record.m * segment.model_record.sigma.powi(3) * count; @@ -252,28 +225,20 @@ impl ParameterHetero for GcPcSaftEosParameters { }); // Association - let sigma3_kappa_aibj = Array2::from_shape_fn([kappa_ab.len(); 2], |(i, j)| { - (sigma[assoc_segment[i]] * sigma[assoc_segment[j]]).powf(1.5) - * (kappa_ab[i] * kappa_ab[j]).sqrt() - }); - let epsilon_k_aibj = Array2::from_shape_fn([epsilon_k_ab.len(); 2], |(i, j)| { - 0.5 * (epsilon_k_ab[i] + epsilon_k_ab[j]) - }); + let sigma = Array1::from_vec(sigma); + let component_index = Array1::from_vec(component_index); + let association = + AssociationParameters::new(&association_records, &sigma, Some(&component_index)); Ok(Self { molarweight, - component_index: Array1::from_vec(component_index), + component_index, identifiers, - n: Array1::from_vec(n), m: Array1::from_vec(m), - sigma: Array1::from_vec(sigma), + sigma, epsilon_k: Array1::from_vec(epsilon_k), bonds, - assoc_segment: Array1::from_vec(assoc_segment), - kappa_ab: Array1::from_vec(kappa_ab), - epsilon_k_ab: Array1::from_vec(epsilon_k_ab), - na: Array1::from_vec(na), - nb: Array1::from_vec(nb), + association, dipole_comp: Array1::from_vec(dipole_comp), mu: Array1::from_vec(mu), mu2: Array1::from_vec(mu2), @@ -283,8 +248,6 @@ impl ParameterHetero for GcPcSaftEosParameters { k_ij, sigma_ij, epsilon_k_ij, - sigma3_kappa_aibj, - epsilon_k_aibj, chemical_records, segment_records, binary_segment_records, @@ -327,7 +290,7 @@ impl GcPcSaftEosParameters { let o = &mut output; write!( o, - "|component|molarweight|dipole moment|segment|count|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|$\\mu$|$Q$|\n|-|-|-|-|-|-|-|-|-|-|-|-|-|-|" + "|component|molarweight|dipole moment|group|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); for i in 0..self.m.len() { @@ -339,7 +302,7 @@ impl GcPcSaftEosParameters { "{}|{}|{}", pure.name .as_ref() - .unwrap_or(&format!("Component {}", i + 1)), + .unwrap_or(&format!("Component {}", self.component_index[i] + 1)), self.molarweight[self.component_index[i]], if let Some(d) = self.dipole_comp.iter().position(|&d| d == i) { format!("{}", self.mu[d]) @@ -348,20 +311,23 @@ impl GcPcSaftEosParameters { } ) }; - let association = if let Some(a) = self.assoc_segment.iter().position(|&a| a == i) { - format!( - "{}|{}|{}|{}", - self.kappa_ab[a], self.epsilon_k_ab[a], self.na[a], self.nb[a] - ) - } else { - "|||".to_string() - }; + let association = + if let Some(a) = self.association.assoc_comp.iter().position(|&a| a == i) { + format!( + "{}|{}|{}|{}", + self.association.kappa_ab[a], + self.association.epsilon_k_ab[a], + self.association.na[a], + self.association.nb[a] + ) + } else { + "|||".to_string() + }; write!( o, - "\n|{}|{}|{}|{}|{}|{}|{}|||", + "\n|{}|{}|{}|{}|{}|{}|", component, self.identifiers[i], - self.n[i], self.m[i], self.sigma[i], self.epsilon_k[i], @@ -369,7 +335,7 @@ impl GcPcSaftEosParameters { ) .unwrap(); } - write!(o, "\n\n|component|segment 1|segment 2|bonds|\n|-|-|-|-|").unwrap(); + write!(o, "\n\n|component|group 1|group 2|bonds|\n|-|-|-|-|").unwrap(); let mut last_component = None; for ([c1, c2], &c) in &self.bonds { @@ -411,13 +377,13 @@ impl std::fmt::Display for GcPcSaftEosParameters { write!(f, "\n\tsigma={}", self.sigma)?; write!(f, "\n\tepsilon_k={}", self.epsilon_k)?; write!(f, "\n\tbonds={:?}", self.bonds)?; - if !self.assoc_segment.is_empty() { - write!(f, "\n\tassoc_segment={}", self.assoc_segment)?; - write!(f, "\n\tkappa_ab={}", self.kappa_ab)?; - write!(f, "\n\tepsilon_k_ab={}", self.epsilon_k_ab)?; - write!(f, "\n\tna={}", self.na)?; - write!(f, "\n\tnb={}", self.nb)?; - } + // if !self.assoc_segment.is_empty() { + // write!(f, "\n\tassoc_segment={}", self.assoc_segment)?; + // write!(f, "\n\tkappa_ab={}", self.kappa_ab)?; + // write!(f, "\n\tepsilon_k_ab={}", self.epsilon_k_ab)?; + // write!(f, "\n\tna={}", self.na)?; + // write!(f, "\n\tnb={}", self.nb)?; + // } if !self.dipole_comp.is_empty() { write!(f, "\n\tdipole_comp={}", self.dipole_comp)?; write!(f, "\n\tmu={}", self.mu)?; @@ -430,12 +396,13 @@ impl std::fmt::Display for GcPcSaftEosParameters { pub mod test { use super::*; use feos_core::parameter::{ChemicalRecord, Identifier}; + use feos_saft::AssociationRecord; fn ch3() -> SegmentRecord { SegmentRecord::new( "CH3".into(), 15.0, - GcPcSaftRecord::new(0.77247, 3.6937, 181.49, None, None, None, None, None, None), + GcPcSaftRecord::new(0.77247, 3.6937, 181.49, None, None, None), None, ) } @@ -444,7 +411,7 @@ pub mod test { SegmentRecord::new( "CH2".into(), 14.0, - GcPcSaftRecord::new(0.7912, 3.0207, 157.23, None, None, None, None, None, None), + GcPcSaftRecord::new(0.7912, 3.0207, 157.23, None, None, None), None, ) } @@ -458,10 +425,7 @@ pub mod test { 2.7702, 334.29, None, - Some(0.009583), - Some(2575.9), - None, - None, + Some(AssociationRecord::new(0.009583, 2575.9, 1.0, 1.0)), None, ), None, diff --git a/src/gc_pcsaft/mod.rs b/src/gc_pcsaft/mod.rs index 96f67262c..b42e54a4d 100644 --- a/src/gc_pcsaft/mod.rs +++ b/src/gc_pcsaft/mod.rs @@ -3,7 +3,7 @@ #[cfg(feature = "dft")] mod dft; -mod eos; +pub(crate) mod eos; #[cfg(feature = "micelles")] pub mod micelles; mod record; diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs index 5e25c9ca7..3b874a6a2 100644 --- a/src/gc_pcsaft/python/mod.rs +++ b/src/gc_pcsaft/python/mod.rs @@ -2,6 +2,7 @@ use super::dft::GcPcSaftFunctionalParameters; use super::eos::GcPcSaftEosParameters; use super::record::GcPcSaftRecord; +use crate::association::python::PyAssociationRecord; use feos_core::joback::JobackRecord; use feos_core::parameter::{ BinaryRecord, IdentifierOption, ParameterError, ParameterHetero, SegmentRecord, @@ -18,9 +19,7 @@ use std::rc::Rc; mod micelles; #[pyclass(name = "GcPcSaftRecord", unsendable)] -#[pyo3( - text_signature = "(m, sigma, epsilon_k, mu=None, q=None, kappa_ab=None, epsilon_k_ab=None, na=None, nb=None)" -)] +#[pyo3(text_signature = "(m, sigma, epsilon_k, mu=None, association_record=None, psi_dft=None)")] #[derive(Clone)] pub struct PyGcPcSaftRecord(GcPcSaftRecord); @@ -32,10 +31,7 @@ impl PyGcPcSaftRecord { sigma: f64, epsilon_k: f64, mu: Option, - kappa_ab: Option, - epsilon_k_ab: Option, - na: Option, - nb: Option, + association_record: Option, psi_dft: Option, ) -> Self { Self(GcPcSaftRecord::new( @@ -43,14 +39,36 @@ impl PyGcPcSaftRecord { sigma, epsilon_k, mu, - kappa_ab, - epsilon_k_ab, - na, - nb, + association_record.map(|r| r.0), psi_dft, )) } + #[getter] + fn get_m(&self) -> f64 { + self.0.m + } + + #[getter] + fn get_sigma(&self) -> f64 { + self.0.sigma + } + + #[getter] + fn get_epsilon_k(&self) -> f64 { + self.0.epsilon_k + } + + #[getter] + fn get_mu(&self) -> Option { + self.0.mu + } + + #[getter] + fn get_association_record(&self) -> Option { + self.0.association_record.clone().map(PyAssociationRecord) + } + fn __repr__(&self) -> PyResult { Ok(self.0.to_string()) } @@ -99,9 +117,9 @@ impl_parameter_from_segments!(GcPcSaftFunctionalParameters, PyGcPcSaftFunctional #[cfg(feature = "dft")] #[pymethods] impl PyGcPcSaftFunctionalParameters { - fn _repr_markdown_(&self) -> String { - self.0.to_markdown() - } + // fn _repr_markdown_(&self) -> String { + // self.0.to_markdown() + // } #[getter] fn get_graph(&self, py: Python) -> PyResult { @@ -138,6 +156,7 @@ pub fn gc_pcsaft(_py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + #[cfg(feature = "dft")] m.add_class::()?; Ok(()) } diff --git a/src/gc_pcsaft/record.rs b/src/gc_pcsaft/record.rs index ead9dce86..6f0a849a0 100644 --- a/src/gc_pcsaft/record.rs +++ b/src/gc_pcsaft/record.rs @@ -1,7 +1,8 @@ +use feos_saft::AssociationRecord; use serde::{Deserialize, Serialize}; /// gc-PC-SAFT pure-component parameters. -#[derive(Serialize, Deserialize, Debug, Clone, Default)] +#[derive(Serialize, Deserialize, Clone, Default)] pub struct GcPcSaftRecord { /// Segment shape factor pub m: f64, @@ -10,27 +11,12 @@ pub struct GcPcSaftRecord { /// Energetic parameter in units of Kelvin pub epsilon_k: f64, /// Dipole moment in units of Debye - #[serde(default)] #[serde(skip_serializing_if = "Option::is_none")] pub mu: Option, - /// association volume parameter - #[serde(default)] + /// Association parameters #[serde(skip_serializing_if = "Option::is_none")] - pub kappa_ab: Option, - /// association energy parameter - #[serde(default)] - #[serde(skip_serializing_if = "Option::is_none")] - pub epsilon_k_ab: Option, - /// \# of association sites of type A - #[serde(default)] - #[serde(skip_serializing_if = "Option::is_none")] - pub na: Option, - /// \# of association sites of type B - #[serde(default)] - #[serde(skip_serializing_if = "Option::is_none")] - pub nb: Option, + pub association_record: Option, /// interaction range parameter for the dispersion functional - #[serde(default)] #[serde(skip_serializing_if = "Option::is_none")] pub psi_dft: Option, } @@ -41,10 +27,7 @@ impl GcPcSaftRecord { sigma: f64, epsilon_k: f64, mu: Option, - kappa_ab: Option, - epsilon_k_ab: Option, - na: Option, - nb: Option, + association_record: Option, psi_dft: Option, ) -> Self { Self { @@ -52,10 +35,7 @@ impl GcPcSaftRecord { sigma, epsilon_k, mu, - kappa_ab, - epsilon_k_ab, - na, - nb, + association_record, psi_dft, } } @@ -69,17 +49,8 @@ impl std::fmt::Display for GcPcSaftRecord { if let Some(n) = &self.mu { write!(f, ", mu={}", n)?; } - if let Some(n) = &self.kappa_ab { - write!(f, ", kappa_ab={}", n)?; - } - if let Some(n) = &self.epsilon_k_ab { - write!(f, ", epsilon_k_ab={}", n)?; - } - if let Some(n) = &self.na { - write!(f, ", na={}", n)?; - } - if let Some(n) = &self.nb { - write!(f, ", nb={}", n)?; + if let Some(n) = &self.association_record { + write!(f, ", association_record={}", n)?; } write!(f, ")") } diff --git a/src/pcsaft/dft/association.rs b/src/pcsaft/dft/association.rs deleted file mode 100644 index a5fadf396..000000000 --- a/src/pcsaft/dft/association.rs +++ /dev/null @@ -1,188 +0,0 @@ -use super::PcSaftParameters; -use crate::pcsaft::eos::association::{ - assoc_site_frac_a, assoc_site_frac_ab, helmholtz_energy_density_cross_association, -}; -use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; -use feos_saft::HardSphereProperties; -use ndarray::*; -use num_dual::DualNum; -use std::f64::consts::PI; -use std::fmt; -use std::rc::Rc; - -pub const N0_CUTOFF: f64 = 1e-9; - -#[derive(Clone)] -pub struct AssociationFunctional { - parameters: Rc, - max_iter: usize, - tol: f64, -} - -impl AssociationFunctional { - pub fn new(parameters: Rc, max_iter: usize, tol: f64) -> Self { - Self { - parameters, - max_iter, - tol, - } - } -} - -impl FunctionalContributionDual for AssociationFunctional -where - N: DualNum + ScalarOperand, -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { - let p = &self.parameters; - let r = p.hs_diameter(temperature) * 0.5; - WeightFunctionInfo::new(Array1::from_shape_fn(r.len(), |i| i), false) - .add( - WeightFunction::new_scaled(r.clone(), WeightFunctionShape::Delta), - false, - ) - .add( - WeightFunction { - prefactor: p.m.mapv(N::from), - kernel_radius: r.clone(), - shape: WeightFunctionShape::DeltaVec, - }, - false, - ) - .add( - WeightFunction { - prefactor: p.m.mapv(N::from), - kernel_radius: r, - shape: WeightFunctionShape::Theta, - }, - true, - ) - } - - fn calculate_helmholtz_energy_density( - &self, - temperature: N, - weighted_densities: ArrayView2, - ) -> Result, EosError> { - let p = &self.parameters; - // number of components - let n = p.m.len(); - // number of dimensions - let dim = (weighted_densities.shape()[0] - 1) / n - 1; - - // weighted densities - let n0i = weighted_densities.slice_axis(Axis(0), Slice::new(0, Some(n as isize), 1)); - let n2vi: Vec<_> = (0..dim) - .into_iter() - .map(|i| { - weighted_densities.slice_axis( - Axis(0), - Slice::new((n * (i + 1)) as isize, Some((n * (i + 2)) as isize), 1), - ) - }) - .collect(); - let n3 = weighted_densities.index_axis(Axis(0), n * (dim + 1)); - - // calculate rho0 - let r = p.hs_diameter(temperature) * 0.5; - let mut n2i = Array::zeros(n0i.raw_dim()); - for i in 0..n { - n2i.index_axis_mut(Axis(0), i) - .assign(&(&n0i.index_axis(Axis(0), i) * (r[i].powi(2) * (p.m[i] * 4.0 * PI)))); - } - let mut rho0: Array2 = (n2vi - .iter() - .map(|n2vi| n2vi * n2vi) - .fold(Array::zeros(n0i.raw_dim()), |acc, x| acc + x) - / -(&n2i * &n2i) - + 1.0) - * n0i; - rho0.iter_mut().zip(&n0i).for_each(|(rho0, &n0i)| { - if n0i.re() < N0_CUTOFF { - *rho0 = n0i; - } - }); - - // calculate xi - let n2v: Vec<_> = n2vi.iter().map(|n2vi| n2vi.sum_axis(Axis(0))).collect(); - let n2 = n2i.sum_axis(Axis(0)); - let mut xi = n2v - .iter() - .map(|n2v| n2v * n2v) - .fold(Array::zeros(n3.raw_dim()), |acc, x| acc + x) - / -(&n2 * &n2) - + 1.0; - xi.iter_mut() - .zip(&n0i.sum_axis(Axis(0))) - .for_each(|(xi, &n0i)| { - if n0i.re() < N0_CUTOFF { - *xi = N::one(); - } - }); - - // auxiliary variables - let n3i = n3.mapv(|n3| (-n3 + 1.0).recip()); - - // only one associating component - if p.nassoc == 1 { - // association strength - let a = p.assoc_comp[0]; - let k = &n2 * &n3i * r[a]; - let deltarho = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) - * ((temperature.recip() * p.epsilon_k_aibj[(a, a)]).exp_m1() - * (p.sigma[a].powi(3) * p.kappa_aibj[(a, a)])) - * rho0.index_axis(Axis(0), a); - - let na = p.na[a]; - let nb = p.nb[a]; - let f = |x: N| x.ln() - x * 0.5 + 0.5; - if nb > 0.0 { - // no cross association, two association sites - let xa = deltarho.mapv(|d| assoc_site_frac_ab(d, na, nb)); - let xb = (&xa - 1.0) * (na / nb) + 1.0; - Ok((xa.mapv(f) * na + xb.mapv(f) * nb) * rho0.index_axis(Axis(0), a)) - } else { - // no cross association, one association site - let xa = deltarho.mapv(|d| assoc_site_frac_a(d, na)); - - Ok(xa.mapv(f) * na * rho0.index_axis(Axis(0), a)) - } - } else { - let mut x: Array1 = Array::from_elem(2 * p.nassoc, 0.2); - Ok(rho0 - .view() - .into_shape([n, rho0.len() / n]) - .unwrap() - .axis_iter(Axis(1)) - .zip(n2.iter()) - .zip(n3i.iter()) - .zip(xi.iter()) - .map(|(((rho0, &n2), &n3i), &xi)| { - helmholtz_energy_density_cross_association( - p, - temperature, - &rho0, - &r, - n2, - n3i, - xi, - self.max_iter, - self.tol, - Some(&mut x), - ) - }) - .collect::, _>>()? - .into_shape(n2.raw_dim()) - .unwrap()) - } - } -} - -impl fmt::Display for AssociationFunctional { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Association functional") - } -} diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index 6e949a948..257597dfd 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -1,27 +1,25 @@ use super::PcSaftParameters; use crate::pcsaft::eos::PcSaftOptions; -use association::AssociationFunctional; -use dispersion::AttractiveFunctional; use feos_core::joback::Joback; use feos_core::parameter::Parameter; use feos_core::{IdealGasContribution, MolarWeight}; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use feos_saft::{FMTContribution, FMTVersion}; -use hard_chain::ChainFunctional; +use feos_saft::{Association, FMTContribution, FMTVersion}; use ndarray::{Array1, Array2}; use num_traits::One; -use pure_saft_functional::*; use quantity::si::*; use std::f64::consts::FRAC_PI_6; use std::rc::Rc; -mod association; mod dispersion; mod hard_chain; mod polar; mod pure_saft_functional; +use dispersion::AttractiveFunctional; +use hard_chain::ChainFunctional; +use pure_saft_functional::*; /// PC-SAFT Helmholtz energy functional. pub struct PcSaftFunctional { @@ -77,9 +75,10 @@ impl PcSaftFunctional { contributions.push(Box::new(att)); // Association - if parameters.nassoc > 0 { - let assoc = AssociationFunctional::new( - parameters.clone(), + if !parameters.association.assoc_comp.is_empty() { + let assoc = Association::new( + ¶meters, + ¶meters.association, saft_options.max_iter_cross_assoc, saft_options.tol_cross_assoc, ); diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index d60dd88d0..85098cff1 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -1,14 +1,13 @@ -use super::association::N0_CUTOFF; use super::polar::{pair_integral_ij, triplet_integral_ijk}; use super::PcSaftParameters; -use crate::pcsaft::eos::association::{assoc_site_frac_a, assoc_site_frac_ab}; +// use crate::association::dft::N0_CUTOFF; use crate::pcsaft::eos::dispersion::{A0, A1, A2, B0, B1, B2}; use crate::pcsaft::eos::polar::{AD, AQ, BD, BQ, CD, CQ, PI_SQ_43}; use feos_core::{EosError, EosResult}; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; -use feos_saft::{FMTVersion, HardSphereProperties}; +use feos_saft::{Association, FMTVersion, HardSphereProperties}; use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; @@ -17,6 +16,7 @@ use std::rc::Rc; const PI36M1: f64 = 1.0 / (36.0 * PI); const N3_CUTOFF: f64 = 1e-5; +const N0_CUTOFF: f64 = 1e-9; #[derive(Clone)] pub struct PureFMTAssocFunctional { @@ -112,7 +112,8 @@ impl + ScalarOperand> FunctionalContributionDual for PureFMTA let mut phi = -(&n0 * &ln31) + n1n2 * &n3m1rec + n2n2 * n2 * PI36M1 * f3; // association - if p.nassoc == 1 { + let a = &p.association; + if a.assoc_comp.len() == 1 { let mut xi = -(&n2v * &n2v).sum_axis(Axis(0)) / (&n2 * &n2) + 1.0; xi.iter_mut().zip(&n2).for_each(|(xi, &n2)| { if n2.re() < N0_CUTOFF * 4.0 * PI * p.m[0] * r.re().powi(2) { @@ -122,19 +123,22 @@ impl + ScalarOperand> FunctionalContributionDual for PureFMTA let k = &n2 * &n3m1rec * r; let deltarho = (((&k / 18.0 + 0.5) * &k * &xi + 1.0) * n3m1rec) - * ((temperature.recip() * p.epsilon_k_aibj[(0, 0)]).exp_m1() - * (p.sigma[0].powi(3) * p.kappa_aibj[(0, 0)])) + * ((temperature.recip() * a.epsilon_k_aibj[(0, 0)]).exp_m1() + * a.sigma3_kappa_aibj[(0, 0)]) * (&n0 / p.m[0] * &xi); let f = |x: N| x.ln() - x * 0.5 + 0.5; phi = phi - + if p.nb[0] > 0.0 { - let xa = deltarho.mapv(|d| assoc_site_frac_ab(d, p.na[0], p.nb[0])); - let xb = (xa.clone() - 1.0) * p.na[0] / p.nb[0] + 1.0; - (n0 / p.m[0] * xi) * (xa.mapv(f) * p.na[0] + xb.mapv(f) * p.nb[0]) + + if a.nb[0] > 0.0 { + let xa = deltarho.mapv(|d| { + Association::::assoc_site_frac_ab(d, a.na[0], a.nb[0]) + }); + let xb = (xa.clone() - 1.0) * a.na[0] / a.nb[0] + 1.0; + (n0 / p.m[0] * xi) * (xa.mapv(f) * a.na[0] + xb.mapv(f) * a.nb[0]) } else { - let xa = deltarho.mapv(|d| assoc_site_frac_a(d, p.na[0])); - n0 / p.m[0] * xi * (xa.mapv(f) * p.na[0]) + let xa = deltarho + .mapv(|d| Association::::assoc_site_frac_a(d, a.na[0])); + n0 / p.m[0] * xi * (xa.mapv(f) * a.na[0]) }; } diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index 2af502e7e..d211c6500 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -5,21 +5,20 @@ use feos_core::{ Contributions, EntropyScaling, EosError, EosResult, EquationOfState, HelmholtzEnergy, IdealGasContribution, MolarWeight, State, }; -use feos_saft::HardSphere; +use feos_saft::{Association, HardSphere}; use ndarray::Array1; use quantity::si::*; use std::f64::consts::{FRAC_PI_6, PI}; use std::rc::Rc; -pub(crate) mod association; pub(crate) mod dispersion; pub(crate) mod hard_chain; pub(crate) mod polar; mod qspr; -use association::{Association, CrossAssociation}; use dispersion::Dispersion; use hard_chain::HardChain; -use polar::{DQVariants, Dipole, DipoleQuadrupole, Quadrupole}; +pub use polar::DQVariants; +use polar::{Dipole, DipoleQuadrupole, Quadrupole}; use qspr::QSPR; #[allow(clippy::upper_case_acronyms)] @@ -88,16 +87,13 @@ impl PcSaft { variant: options.dq_variant, })); }; - match parameters.nassoc { - 0 => (), - 1 => contributions.push(Box::new(Association { - parameters: parameters.clone(), - })), - _ => contributions.push(Box::new(CrossAssociation { - parameters: parameters.clone(), - max_iter: options.max_iter_cross_assoc, - tol: options.tol_cross_assoc, - })), + if !parameters.association.assoc_comp.is_empty() { + contributions.push(Box::new(Association::new( + ¶meters, + ¶meters.association, + options.max_iter_cross_assoc, + options.tol_cross_assoc, + ))); }; let joback_records = parameters.joback_records.clone(); diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs index e662ce29e..e0946c7b3 100644 --- a/src/pcsaft/eos/polar.rs +++ b/src/pcsaft/eos/polar.rs @@ -375,7 +375,9 @@ impl fmt::Display for Quadrupole { } } +/// Different combination rules used in the dipole-quadrupole contribution. #[derive(Clone, Copy)] +#[cfg_attr(feature = "python", pyo3::pyclass)] pub enum DQVariants { DQ35, DQ44, diff --git a/src/pcsaft/eos/qspr.rs b/src/pcsaft/eos/qspr.rs index 020133a32..98cf41975 100644 --- a/src/pcsaft/eos/qspr.rs +++ b/src/pcsaft/eos/qspr.rs @@ -70,7 +70,7 @@ pub struct QSPR { impl> IdealGasContributionDual for QSPR { fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1 { - let (c_300, c_400) = match self.parameters.nassoc { + let (c_300, c_400) = match self.parameters.association.assoc_comp.len() { 0 => match self.parameters.ndipole + self.parameters.nquadpole { 0 => (NA_NP_300, NA_NP_400), _ => (NA_P_300, NA_P_400), @@ -85,10 +85,13 @@ impl> IdealGasContributionDual for QSPR { let p1 = epsilon_kt * self.parameters.m[i]; let p2 = sigma3 * self.parameters.m[i]; let p3 = epsilon_kt * p2; - let p4 = (temperature.recip() * self.parameters.epsilon_k_ab[i]).exp_m1() - * p2 - * sigma3 - * self.parameters.kappa_ab[i]; + let p4 = self.parameters.pure_records[i] + .model_record + .association_record + .as_ref() + .map_or(D::zero(), |a| { + (temperature.recip() * a.epsilon_k_ab).exp_m1() * p2 * sigma3 * a.kappa_ab + }); let p5 = p2 * self.parameters.q[i]; let p6 = 1.0; diff --git a/src/pcsaft/mod.rs b/src/pcsaft/mod.rs index bf50c486a..b7415a2f8 100644 --- a/src/pcsaft/mod.rs +++ b/src/pcsaft/mod.rs @@ -4,11 +4,11 @@ #[cfg(feature = "dft")] mod dft; mod eos; -mod parameters; +pub(crate) mod parameters; #[cfg(feature = "dft")] pub use dft::PcSaftFunctional; -pub use eos::{PcSaft, PcSaftOptions}; +pub use eos::{DQVariants, PcSaft, PcSaftOptions}; pub use parameters::{PcSaftParameters, PcSaftRecord}; #[cfg(feature = "python")] diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 3ee93a3c8..b22025233 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -3,7 +3,7 @@ use feos_core::joback::JobackRecord; use feos_core::parameter::{ FromSegments, FromSegmentsBinary, Parameter, ParameterError, PureRecord, }; -use feos_saft::{HardSphereProperties, MonomerShape}; +use feos_saft::{AssociationParameters, AssociationRecord, HardSphereProperties, MonomerShape}; use ndarray::{Array, Array1, Array2}; use num_dual::DualNum; use num_traits::Zero; @@ -27,18 +27,9 @@ pub struct PcSaftRecord { /// Quadrupole moment in units of Debye #[serde(skip_serializing_if = "Option::is_none")] pub q: Option, - /// Association volume parameter + /// Association parameters #[serde(skip_serializing_if = "Option::is_none")] - pub kappa_ab: Option, - /// Association energy parameter in units of Kelvin - #[serde(skip_serializing_if = "Option::is_none")] - pub epsilon_k_ab: Option, - /// \# of association sites of type A - #[serde(skip_serializing_if = "Option::is_none")] - pub na: Option, - /// \# of association sites of type B - #[serde(skip_serializing_if = "Option::is_none")] - pub nb: Option, + pub association_record: Option, /// Entropy scaling coefficients for the viscosity #[serde(skip_serializing_if = "Option::is_none")] pub viscosity: Option<[f64; 4]>, @@ -70,22 +61,22 @@ impl FromSegments for PcSaftRecord { .iter() .filter_map(|(s, n)| s.mu.map(|mu| mu * n)) .reduce(|a, b| a + b); - let kappa_ab = segments - .iter() - .filter_map(|(s, n)| s.kappa_ab.map(|k| k * n)) - .reduce(|a, b| a + b); - let epsilon_k_ab = segments - .iter() - .filter_map(|(s, n)| s.epsilon_k_ab.map(|e| e * n)) - .reduce(|a, b| a + b); - let na = segments - .iter() - .filter_map(|(s, n)| s.na.map(|na| na * n)) - .reduce(|a, b| a + b); - let nb = segments + let association_record = segments .iter() - .filter_map(|(s, n)| s.nb.map(|nb| nb * n)) - .reduce(|a, b| a + b); + .filter_map(|(s, n)| { + s.association_record.as_ref().map(|record| { + [ + record.kappa_ab * n, + record.epsilon_k_ab * n, + record.na * n, + record.nb * n, + ] + }) + }) + .reduce(|a, b| [a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]]) + .map(|[kappa_ab, epsilon_k_ab, na, nb]| { + AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb) + }); // entropy scaling let mut viscosity = if segments @@ -147,10 +138,7 @@ impl FromSegments for PcSaftRecord { epsilon_k: epsilon_k / m, mu, q, - kappa_ab, - epsilon_k_ab, - na, - nb, + association_record, viscosity, diffusion, thermal_conductivity, @@ -173,18 +161,9 @@ impl FromSegments for PcSaftRecord { "{dipole_comps} segment with dipole moment." ))); }; - let faulty_assoc_comps = segments - .iter() - .filter_map(|(s, _)| s.kappa_ab.xor(s.epsilon_k_ab)) - .count(); - if faulty_assoc_comps > 0 { - return Err(ParameterError::IncompatibleParameters(format!( - "Incorrectly specified association sites on {faulty_assoc_comps} segment(s)" - ))); - } let assoc_comps = segments .iter() - .filter_map(|(s, _)| s.kappa_ab.and(s.epsilon_k_ab)) + .filter_map(|(s, _)| s.association_record.as_ref()) .count(); if assoc_comps > 1 { return Err(ParameterError::IncompatibleParameters(format!( @@ -211,17 +190,8 @@ impl std::fmt::Display for PcSaftRecord { if let Some(n) = &self.q { write!(f, ", q={}", n)?; } - if let Some(n) = &self.kappa_ab { - write!(f, ", kappa_ab={}", n)?; - } - if let Some(n) = &self.epsilon_k_ab { - write!(f, ", epsilon_k_ab={}", n)?; - } - if let Some(n) = &self.na { - write!(f, ", na={}", n)?; - } - if let Some(n) = &self.nb { - write!(f, ", nb={}", n)?; + if let Some(n) = &self.association_record { + write!(f, ", association_record={}", n)?; } if let Some(n) = &self.viscosity { write!(f, ", viscosity={:?}", n)?; @@ -243,10 +213,7 @@ impl PcSaftRecord { epsilon_k: f64, mu: Option, q: Option, - kappa_ab: Option, - epsilon_k_ab: Option, - na: Option, - nb: Option, + association_record: Option, viscosity: Option<[f64; 4]>, diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, @@ -257,10 +224,7 @@ impl PcSaftRecord { epsilon_k, mu, q, - kappa_ab, - epsilon_k_ab, - na, - nb, + association_record, viscosity, diffusion, thermal_conductivity, @@ -311,22 +275,15 @@ pub struct PcSaftParameters { pub q: Array1, pub mu2: Array1, pub q2: Array1, - pub kappa_ab: Array1, - pub epsilon_k_ab: Array1, - pub na: Array1, - pub nb: Array1, - pub kappa_aibj: Array2, - pub epsilon_k_aibj: Array2, + pub association: AssociationParameters, pub k_ij: Array2, pub sigma_ij: Array2, pub epsilon_k_ij: Array2, pub e_k_ij: Array2, pub ndipole: usize, pub nquadpole: usize, - pub nassoc: usize, pub dipole_comp: Array1, pub quadpole_comp: Array1, - pub assoc_comp: Array1, pub viscosity: Option>, pub diffusion: Option>, pub thermal_conductivity: Option>, @@ -352,10 +309,7 @@ impl Parameter for PcSaftParameters { let mut epsilon_k = Array::zeros(n); let mut mu = Array::zeros(n); let mut q = Array::zeros(n); - let mut na = Array::zeros(n); - let mut nb = Array::zeros(n); - let mut kappa_ab = Array::zeros(n); - let mut epsilon_k_ab = Array::zeros(n); + let mut association_records = Vec::with_capacity(n); let mut viscosity = Vec::with_capacity(n); let mut diffusion = Vec::with_capacity(n); let mut thermal_conductivity = Vec::with_capacity(n); @@ -370,10 +324,7 @@ impl Parameter for PcSaftParameters { epsilon_k[i] = r.epsilon_k; mu[i] = r.mu.unwrap_or(0.0); q[i] = r.q.unwrap_or(0.0); - na[i] = r.na.unwrap_or(1.0); - nb[i] = r.nb.unwrap_or(1.0); - kappa_ab[i] = r.kappa_ab.unwrap_or(0.0); - epsilon_k_ab[i] = r.epsilon_k_ab.unwrap_or(0.0); + association_records.push(r.association_record.clone()); viscosity.push(r.viscosity); diffusion.push(r.diffusion); thermal_conductivity.push(r.thermal_conductivity); @@ -398,24 +349,8 @@ impl Parameter for PcSaftParameters { .filter_map(|(i, &q2)| (q2.abs() > 0.0).then(|| i)) .collect(); let nquadpole = quadpole_comp.len(); - let assoc_comp: Array1 = kappa_ab - .iter() - .enumerate() - .filter_map(|(i, &k)| (k.abs() > 0.0).then(|| i)) - .collect(); - let nassoc = assoc_comp.len(); - let mut kappa_aibj = Array::zeros([n, n]); - let mut epsilon_k_aibj = Array::zeros([n, n]); - for i in 0..nassoc { - for j in 0..nassoc { - let ai = assoc_comp[i]; - let bj = assoc_comp[j]; - kappa_aibj[[ai, bj]] = (kappa_ab[ai] * kappa_ab[bj]).sqrt() - * (2.0 * (sigma[ai] * sigma[bj]).sqrt() / (sigma[ai] + sigma[bj])).powf(3.0); - epsilon_k_aibj[[ai, bj]] = 0.5 * (epsilon_k_ab[ai] + epsilon_k_ab[bj]); - } - } + let association = AssociationParameters::new(&association_records, &sigma, None); let k_ij = binary_records.map(|br| br.k_ij); let mut epsilon_k_ij = Array::zeros((n, n)); @@ -474,22 +409,15 @@ impl Parameter for PcSaftParameters { q, mu2, q2, - kappa_ab, - epsilon_k_ab, - na, - nb, - kappa_aibj, - epsilon_k_aibj, + association, k_ij, sigma_ij, epsilon_k_ij, e_k_ij, ndipole, nquadpole, - nassoc, dipole_comp, quadpole_comp, - assoc_comp, viscosity: viscosity_coefficients, diffusion: diffusion_coefficients, thermal_conductivity: thermal_conductivity_coefficients, @@ -531,23 +459,28 @@ impl PcSaftParameters { "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$\\mu$|$Q$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); - for i in 0..self.m.len() { - let component = self.pure_records[i].identifier.name.clone(); + for (i, record) in self.pure_records.iter().enumerate() { + let component = record.identifier.name.clone(); let component = component.unwrap_or(format!("Component {}", i + 1)); + let association = record + .model_record + .association_record + .clone() + .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, 0.0, 0.0)); write!( o, "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", component, - self.molarweight[i], - self.m[i], - self.sigma[i], - self.epsilon_k[i], - self.mu[i], - self.q[i], - self.kappa_ab[i], - self.epsilon_k_ab[i], - self.na[i], - self.nb[i] + record.molarweight, + record.model_record.m, + record.model_record.sigma, + record.model_record.epsilon_k, + record.model_record.mu.unwrap_or(0.0), + record.model_record.q.unwrap_or(0.0), + association.kappa_ab, + association.epsilon_k_ab, + association.na, + association.nb ) .unwrap(); } @@ -569,12 +502,12 @@ impl std::fmt::Display for PcSaftParameters { if !self.quadpole_comp.is_empty() { write!(f, "\n\tq={}", self.q)?; } - if !self.assoc_comp.is_empty() { - write!(f, "\n\tkappa_ab={}", self.kappa_ab)?; - write!(f, "\n\tepsilon_k_ab={}", self.epsilon_k_ab)?; - write!(f, "\n\tna={}", self.na)?; - write!(f, "\n\tnb={}", self.nb)?; - } + // if !self.assoc_comp.is_empty() { + // write!(f, "\n\tkappa_ab={}", self.kappa_ab)?; + // write!(f, "\n\tepsilon_k_ab={}", self.epsilon_k_ab)?; + // write!(f, "\n\tna={}", self.na)?; + // write!(f, "\n\tnb={}", self.nb)?; + // } if !self.k_ij.iter().all(|k| k.is_zero()) { write!(f, "\n\tk_ij=\n{}", self.k_ij)?; } @@ -909,8 +842,12 @@ pub mod utils { "m": 1.065587, "sigma": 3.000683, "epsilon_k": 366.5121, - "kappa_ab": 0.034867983, - "epsilon_k_ab": 2500.6706 + "association_record": { + "kappa_ab": 0.034867983, + "epsilon_k_ab": 2500.6706, + "na": 1.0, + "nb": 1.0 + } }, "molarweight": 18.0152 }"#; diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index 0e749cc33..6f70a484d 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -1,5 +1,5 @@ -use super::eos::polar::DQVariants; use super::parameters::{PcSaftBinaryRecord, PcSaftParameters, PcSaftRecord}; +use crate::association::python::PyAssociationRecord; use feos_core::joback::JobackRecord; use feos_core::parameter::{ BinaryRecord, Identifier, IdentifierOption, Parameter, ParameterError, PureRecord, @@ -14,20 +14,10 @@ use pyo3::prelude::*; use std::convert::{TryFrom, TryInto}; use std::rc::Rc; -impl From<&str> for DQVariants { - fn from(str: &str) -> Self { - match str { - "dq35" => Self::DQ35, - "dq44" => Self::DQ44, - _ => panic!("dq_variant must be either \"dq35\" or \"dq44\""), - } - } -} - /// Create a set of PC-Saft parameters from records. #[pyclass(name = "PcSaftRecord", unsendable)] #[pyo3( - text_signature = "(m, sigma, epsilon_k, mu=None, q=None, kappa_ab=None, epsilon_k_ab=None, na=None, nb=None, viscosity=None, diffusion=None, thermal_conductivity=None)" + text_signature = "(m, sigma, epsilon_k, mu=None, q=None, association_record=None, viscosity=None, diffusion=None, thermal_conductivity=None)" )] #[derive(Clone)] pub struct PyPcSaftRecord(PcSaftRecord); @@ -41,10 +31,7 @@ impl PyPcSaftRecord { epsilon_k: f64, mu: Option, q: Option, - kappa_ab: Option, - epsilon_k_ab: Option, - na: Option, - nb: Option, + association_record: Option, viscosity: Option<[f64; 4]>, diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, @@ -55,10 +42,7 @@ impl PyPcSaftRecord { epsilon_k, mu, q, - kappa_ab, - epsilon_k_ab, - na, - nb, + association_record.map(|r| r.0), viscosity, diffusion, thermal_conductivity, @@ -91,23 +75,8 @@ impl PyPcSaftRecord { } #[getter] - fn get_kappa_ab(&self) -> Option { - self.0.kappa_ab - } - - #[getter] - fn get_epsilon_k_ab(&self) -> Option { - self.0.epsilon_k_ab - } - - #[getter] - fn get_na(&self) -> Option { - self.0.na - } - - #[getter] - fn get_nb(&self) -> Option { - self.0.nb + fn get_association_record(&self) -> Option { + self.0.association_record.clone().map(PyAssociationRecord) } #[getter] diff --git a/src/python/dft.rs b/src/python/dft.rs index f6665d3c2..aae0f7493 100644 --- a/src/python/dft.rs +++ b/src/python/dft.rs @@ -9,7 +9,7 @@ use crate::impl_estimator; #[cfg(feature = "pcsaft")] use crate::pcsaft::python::PyPcSaftParameters; #[cfg(feature = "pcsaft")] -use crate::pcsaft::{PcSaftFunctional, PcSaftOptions}; +use crate::pcsaft::{DQVariants, PcSaftFunctional, PcSaftOptions}; #[cfg(feature = "pets")] use crate::pets::python::PyPetsParameters; #[cfg(feature = "pets")] @@ -218,8 +218,8 @@ impl PyFunctionalVariant { /// Maximum number of iterations for cross association. Defaults to 50. /// tol_cross_assoc : float /// Tolerance for convergence of cross association. Defaults to 1e-10. - /// dq_variant : {'dq35', 'dq44'}, optional - /// Combination rule used in the dipole/quadrupole term. Defaults to 'dq35' + /// dq_variant : DQVariants, optional + /// Combination rule used in the dipole/quadrupole term. Defaults to 'DQVariants.DQ35' /// /// Returns /// ------- @@ -230,7 +230,7 @@ impl PyFunctionalVariant { max_eta = "0.5", max_iter_cross_assoc = "50", tol_cross_assoc = "1e-10", - dq_variant = "\"dq35\"" + dq_variant = "DQVariants::DQ35" )] #[staticmethod] #[pyo3( @@ -242,13 +242,13 @@ impl PyFunctionalVariant { max_eta: f64, max_iter_cross_assoc: usize, tol_cross_assoc: f64, - dq_variant: &str, + dq_variant: DQVariants, ) -> Self { let options = PcSaftOptions { max_eta, max_iter_cross_assoc, tol_cross_assoc, - dq_variant: dq_variant.into(), + dq_variant, }; Self(Rc::new( PcSaftFunctional::with_options(parameters.0, fmt_version, options).into(), diff --git a/src/python/eos.rs b/src/python/eos.rs index 45b5556ed..1ed767120 100644 --- a/src/python/eos.rs +++ b/src/python/eos.rs @@ -11,7 +11,7 @@ use crate::impl_estimator_entropy_scaling; #[cfg(feature = "pcsaft")] use crate::pcsaft::python::PyPcSaftParameters; #[cfg(feature = "pcsaft")] -use crate::pcsaft::{PcSaft, PcSaftOptions}; +use crate::pcsaft::{DQVariants, PcSaft, PcSaftOptions}; #[cfg(feature = "pets")] use crate::pets::python::PyPetsParameters; #[cfg(feature = "pets")] @@ -229,8 +229,8 @@ impl PyEosVariant { /// Maximum number of iterations for cross association. Defaults to 50. /// tol_cross_assoc : float /// Tolerance for convergence of cross association. Defaults to 1e-10. - /// dq_variant : {'dq35', 'dq44'}, optional - /// Combination rule used in the dipole/quadrupole term. Defaults to 'dq35' + /// dq_variant : DQVariants, optional + /// Combination rule used in the dipole/quadrupole term. Defaults to 'DQVariants.DQ35' /// /// Returns /// ------- @@ -242,7 +242,7 @@ impl PyEosVariant { max_eta = "0.5", max_iter_cross_assoc = "50", tol_cross_assoc = "1e-10", - dq_variant = "\"dq35\"" + dq_variant = "DQVariants::DQ35" )] #[staticmethod] #[pyo3( @@ -253,13 +253,13 @@ impl PyEosVariant { max_eta: f64, max_iter_cross_assoc: usize, tol_cross_assoc: f64, - dq_variant: &str, + dq_variant: DQVariants, ) -> Self { let options = PcSaftOptions { max_eta, max_iter_cross_assoc, tol_cross_assoc, - dq_variant: dq_variant.into(), + dq_variant, }; Self(Rc::new(EosVariant::PcSaft(PcSaft::with_options( parameters.0, diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index bc7b318c8..69a865dd2 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -3,8 +3,11 @@ use approx::assert_relative_eq; use feos::gc_pcsaft::{ GcPcSaft, GcPcSaftEosParameters, GcPcSaftFunctional, GcPcSaftFunctionalParameters, + GcPcSaftRecord, +}; +use feos_core::parameter::{ + self, ChemicalRecord, Identifier, IdentifierOption, ParameterHetero, SegmentRecord, }; -use feos_core::parameter::{IdentifierOption, ParameterHetero}; use feos_core::{PhaseEquilibrium, State, StateBuilder, Verbosity}; use feos_dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification}; use feos_dft::interface::PlanarInterface; @@ -105,6 +108,123 @@ fn test_bulk_implementation() -> Result<(), Box> { Ok(()) } +#[test] +fn test_bulk_association() -> Result<(), Box> { + let segment_records = SegmentRecord::from_json("parameters/pcsaft/sauer2014_hetero.json")?; + let ethylene_glycol = ChemicalRecord::new( + Identifier::default(), + vec!["OH".into(), "CH2".into(), "CH2".into(), "OH".into()], + None, + ); + let eos_parameters = Rc::new(GcPcSaftEosParameters::from_segments( + vec![ethylene_glycol.clone()], + segment_records.clone(), + None, + )?); + let eos = Rc::new(GcPcSaft::new(eos_parameters.clone())); + let func_parameters = Rc::new(GcPcSaftFunctionalParameters::from_segments( + vec![ethylene_glycol], + segment_records, + None, + )?); + let func = Rc::new(GcPcSaftFunctional::new(func_parameters.clone())); + + let t = 200.0 * KELVIN; + let v = 0.002 * METER.powi(3); + let n = arr1(&[1.5]) * MOL; + let state_eos = State::new_nvt(&eos, t, v, &n)?; + let state_func = State::new_nvt(&func, t, v, &n)?; + let p_eos = state_eos.pressure_contributions(); + let p_func = state_func.pressure_contributions(); + println!( + "Equation of state: + \tcomps: {} + \tkappa_ab: {} + \tepsilon_k_ab: {} + \tna: {} + \tnb: {}", + eos_parameters.association.assoc_comp, + eos_parameters.association.kappa_ab, + eos_parameters.association.epsilon_k_ab, + eos_parameters.association.na, + eos_parameters.association.nb, + ); + for (s, x) in &p_eos { + println!("{s:18}: {x:21.16}"); + } + println!( + "\nHelmholtz energy functional: + \tcomps: {} + \tkappa_ab: {} + \tepsilon_k_ab: {} + \tna: {} + \tnb: {}", + func_parameters.association.assoc_comp, + func_parameters.association.kappa_ab, + func_parameters.association.epsilon_k_ab, + func_parameters.association.na, + func_parameters.association.nb, + ); + for (s, x) in &p_func { + println!("{s:26}: {x:21.16}"); + } + assert_relative_eq!(p_eos[4].1, p_func[4].1, max_relative = 1e-14); + // println!("{:29}: {}", p_eos[0].0, p_eos[0].1); + // println!("{:29}: {}", p_func[0].0, p_func[0].1); + // println!(); + // println!("{:29}: {}", p_eos[1].0, p_eos[1].1); + // println!("{:29}: {}", p_func[1].0, p_func[1].1); + // println!(); + // println!("{:29}: {}", p_eos[2].0, p_eos[2].1); + // println!("{:29}: {}", p_func[2].0, p_func[2].1); + // println!(); + // println!("{:29}: {}", p_eos[3].0, p_eos[3].1); + // println!("{:29}: {}", p_func[3].0, p_func[3].1); + + // assert_relative_eq!( + // p_eos[0].1, + // 1.2471689792172869 * MEGA * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + // assert_relative_eq!( + // p_eos[1].1, + // 280.0635060891395938 * KILO * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + // assert_relative_eq!( + // p_eos[2].1, + // -141.9023918353318550 * KILO * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + // assert_relative_eq!( + // p_eos[3].1, + // -763.2289230004602132 * KILO * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + + // assert_relative_eq!( + // p_func[0].1, + // 1.2471689792172869 * MEGA * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + // assert_relative_eq!( + // p_func[1].1, + // 280.0635060891395938 * KILO * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + // assert_relative_eq!( + // p_func[2].1, + // -141.9023918353318550 * KILO * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + // assert_relative_eq!( + // p_func[3].1, + // -763.2289230004602132 * KILO * PASCAL * KB / KB_old, + // max_relative = 1e-14, + // ); + Ok(()) +} + #[test] #[allow(non_snake_case)] fn test_dft() -> Result<(), Box> { diff --git a/tests/pcsaft/test_parameters.json b/tests/pcsaft/test_parameters.json index 988fcb3b8..ef4ad70ca 100644 --- a/tests/pcsaft/test_parameters.json +++ b/tests/pcsaft/test_parameters.json @@ -92,8 +92,12 @@ "m": 1.065587, "sigma": 3.000683, "epsilon_k": 366.5121, - "kappa_ab": 0.034867983, - "epsilon_k_ab": 2500.6706 + "association_record": { + "kappa_ab": 0.034867983, + "epsilon_k_ab": 2500.6706, + "na": 1.0, + "nb": 1.0 + } }, "molarweight": 18.0152 }, From b8c31b38da680d15372524e4f5a02d70328eacd0 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sun, 3 Jul 2022 11:58:59 +0200 Subject: [PATCH 3/9] cleanup --- src/gc_pcsaft/dft/mod.rs | 2 -- src/gc_pcsaft/eos/mod.rs | 2 -- 2 files changed, 4 deletions(-) diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 65c999097..10ec4dce0 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -11,11 +11,9 @@ use quantity::si::{SIArray1, SIUnit, GRAM, MOL}; use std::f64::consts::FRAC_PI_6; use std::rc::Rc; -// mod association; mod dispersion; mod hard_chain; mod parameter; -// use association::AssociationFunctional; use dispersion::AttractiveFunctional; use hard_chain::ChainFunctional; pub use parameter::GcPcSaftFunctionalParameters; diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index 09cb0fcd2..edce4c102 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -7,12 +7,10 @@ use quantity::si::*; use std::f64::consts::FRAC_PI_6; use std::rc::Rc; -// pub(crate) mod association; pub(crate) mod dispersion; mod hard_chain; pub(crate) mod parameter; mod polar; -// use association::{Association, CrossAssociation}; use dispersion::Dispersion; use hard_chain::HardChain; pub use parameter::{GcPcSaftChemicalRecord, GcPcSaftEosParameters}; From 234203f88ff5f8db3dbcaff8197f0ba1b1806025 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 4 Jul 2022 09:48:41 +0200 Subject: [PATCH 4/9] remove obsolete files --- src/gc_pcsaft/eos/association.rs | 375 ------------------------------- src/gc_pcsaft/eos/hard_sphere.rs | 129 ----------- src/pcsaft/eos/association.rs | 295 ------------------------ src/pcsaft/eos/hard_sphere.rs | 112 --------- 4 files changed, 911 deletions(-) delete mode 100644 src/gc_pcsaft/eos/association.rs delete mode 100644 src/gc_pcsaft/eos/hard_sphere.rs delete mode 100644 src/pcsaft/eos/association.rs delete mode 100644 src/pcsaft/eos/hard_sphere.rs diff --git a/src/gc_pcsaft/eos/association.rs b/src/gc_pcsaft/eos/association.rs deleted file mode 100644 index 9f9ef169a..000000000 --- a/src/gc_pcsaft/eos/association.rs +++ /dev/null @@ -1,375 +0,0 @@ -use super::GcPcSaftEosParameters; -use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; -use feos_saft::HardSphereProperties; -use ndarray::*; -use num_dual::linalg::{norm, LU}; -use num_dual::*; -use std::fmt; -use std::rc::Rc; - -#[derive(Clone)] -pub struct Association { - pub parameters: Rc, -} - -#[derive(Clone)] -pub struct CrossAssociation { - pub parameters: Rc, - pub max_iter: usize, - pub tol: f64, -} - -fn association_strength>( - assoc_segment: &Array1, - sigma3_kappa_aibj: &Array2, - epsilon_k_aibj: &Array2, - temperature: D, - diameter: &Array1, - n2: D, - n3i: D, - xi: D, - i: usize, - j: usize, -) -> D { - let ai = assoc_segment[i]; - let aj = assoc_segment[j]; - let k = diameter[ai] * diameter[aj] / (diameter[ai] + diameter[aj]) * (n2 * n3i); - n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) - * sigma3_kappa_aibj[(i, j)] - * (temperature.recip() * epsilon_k_aibj[(i, j)]).exp_m1() -} - -impl> HelmholtzEnergyDual for Association { - fn helmholtz_energy(&self, state: &StateHD) -> D { - let p = &self.parameters; - let c = p.component_index[p.assoc_segment[0]]; - - // temperature dependent segment diameter - let diameter = p.hs_diameter(state.temperature); - - // Packing fractions - let [zeta2, n3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); - let n2 = zeta2 * 6.0; - let n3i = (-n3 + 1.0).recip(); - - // association strength - let deltarho = association_strength( - &p.assoc_segment, - &p.sigma3_kappa_aibj, - &p.epsilon_k_aibj, - state.temperature, - &diameter, - n2, - n3i, - D::one(), - 0, - 0, - ) * state.partial_density[c]; - - let na = p.na[0]; - let nb = p.nb[0]; - if nb > 0.0 { - // no cross association, two association sites - let xa = assoc_site_frac_ab(deltarho, na, nb); - let xb = (xa - 1.0) * (na / nb) + 1.0; - state.moles[c] - * p.n[0] - * ((xa.ln() - xa * 0.5 + 0.5) * na + (xb.ln() - xb * 0.5 + 0.5) * nb) - } else { - // no cross association, one association site - let xa = assoc_site_frac_a(deltarho, na); - - state.moles[c] * p.n[0] * (xa.ln() - xa * 0.5 + 0.5) * na - } - } -} - -impl fmt::Display for Association { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Association") - } -} - -pub(crate) fn assoc_site_frac_ab>(deltarho: D, na: f64, nb: f64) -> D { - if deltarho.re() > f64::EPSILON.sqrt() { - (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt() - - (deltarho * (nb - na) + 1.0)) - / (deltarho * na * 2.0) - } else { - D::one() + deltarho * nb * (deltarho * (nb + na) - 1.0) - } -} - -pub(crate) fn assoc_site_frac_a>(deltarho: D, na: f64) -> D { - if deltarho.re() > f64::EPSILON.sqrt() { - ((deltarho * na * 4.0 + 1.0).powi(2) - 1.0).sqrt() / (deltarho * na * 2.0) - } else { - D::one() + deltarho * na * (deltarho * na * 2.0 - 1.0) - } -} - -impl + ScalarOperand> HelmholtzEnergyDual for CrossAssociation { - fn helmholtz_energy(&self, state: &StateHD) -> D { - let p = &self.parameters; - - // temperature dependent segment diameter - let diameter = p.hs_diameter(state.temperature); - - // Packing fractions - let [zeta2, n3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); - let n2 = zeta2 * 6.0; - let n3i = (-n3 + 1.0).recip(); - - // extract densities of associating segments - let rho_assoc = p - .assoc_segment - .mapv(|a| state.partial_density[p.component_index[a]]) - * &p.n; - - helmholtz_energy_density_cross_association( - &p.assoc_segment, - &p.sigma3_kappa_aibj, - &p.epsilon_k_aibj, - &p.na, - &p.nb, - state.temperature, - &rho_assoc, - &diameter, - n2, - n3i, - D::one(), - self.max_iter, - self.tol, - None, - ) - .unwrap_or_else(|_| D::from(std::f64::NAN)) - * state.volume - } -} - -impl fmt::Display for CrossAssociation { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Cross-association") - } -} - -pub fn helmholtz_energy_density_cross_association + ScalarOperand>( - // p: &GcPcSaftParameters, - assoc_segment: &Array1, - sigma3_kappa_aibj: &Array2, - epsilon_k_aibj: &Array2, - na: &Array1, - nb: &Array1, - temperature: D, - density: &ArrayBase, - diameter: &Array1, - n2: D, - n3i: D, - xi: D, - max_iter: usize, - tol: f64, - x0: Option<&mut Array1>, -) -> Result -where - S: Data, -{ - // check if density is close to 0 - if density.sum().re() < f64::EPSILON { - if let Some(x0) = x0 { - x0.fill(1.0); - } - return Ok(D::zero()); - } - - // number of associating segments - let nassoc = assoc_segment.len(); - - // association strength - let delta = Array::from_shape_fn([nassoc; 2], |(i, j)| { - association_strength( - assoc_segment, - sigma3_kappa_aibj, - epsilon_k_aibj, - temperature, - diameter, - n2, - n3i, - xi, - i, - j, - ) - }); - - // cross-association according to Michelsen2006 - // initialize monomer fraction - let mut x = match &x0 { - Some(x0) => (*x0).clone(), - None => Array::from_elem(2 * nassoc, 0.2), - }; - - for k in 0..max_iter { - if newton_step_cross_association::<_, f64>( - &mut x, - nassoc, - &delta.map(D::re), - na, - nb, - &density.map(D::re), - tol, - )? { - break; - } - if k == max_iter - 1 { - return Err(EosError::NotConverged("Cross association".into())); - } - } - - // calculate derivatives - let mut x_dual = x.mapv(D::from); - for _ in 0..D::NDERIV { - newton_step_cross_association(&mut x_dual, nassoc, &delta, na, nb, density, tol)?; - } - - // save monomer fraction - if let Some(x0) = x0 { - *x0 = x; - } - - // Helmholtz energy density - let xa = x_dual.slice(s![..nassoc]); - let xb = x_dual.slice(s![nassoc..]); - let f = |x: D| x.ln() - x * 0.5 + 0.5; - Ok((density * (xa.mapv(f) * na + xb.mapv(f) * nb)).sum()) -} - -fn newton_step_cross_association + ScalarOperand>( - x: &mut Array1, - nassoc: usize, - delta: &Array2, - na: &Array1, - nb: &Array1, - rho: &ArrayBase, - tol: f64, -) -> Result -where - S: Data, -{ - // gradient - let mut g: Array1 = Array::zeros(2 * nassoc); - // Hessian - let mut h: Array2 = Array::zeros((2 * nassoc, 2 * nassoc)); - - // slice arrays - let (xa, xb) = x.multi_slice_mut((s![..nassoc], s![nassoc..])); - let (mut ga, mut gb) = g.multi_slice_mut((s![..nassoc], s![nassoc..])); - let (mut haa, mut hab, mut hba, mut hbb) = h.multi_slice_mut(( - s![..nassoc, ..nassoc], - s![..nassoc, nassoc..], - s![nassoc.., ..nassoc], - s![nassoc.., nassoc..], - )); - - // calculate gradients and approximate Hessian - for i in 0..nassoc { - let d = &delta.index_axis(Axis(0), i) * rho; - - let dnx = (&xb * nb * &d).sum() + 1.0; - ga[i] = xa[i].recip() - dnx; - hab.index_axis_mut(Axis(0), i).assign(&(&d * &(-nb))); - haa[(i, i)] = -dnx / xa[i]; - - let dnx = (&xa * na * &d).sum() + 1.0; - gb[i] = xb[i].recip() - dnx; - hba.index_axis_mut(Axis(0), i).assign(&(&d * &(-na))); - hbb[(i, i)] = -dnx / xb[i]; - } - - // Newton step - x.assign(&(&*x - &LU::new(h)?.solve(&g))); - - // check convergence - Ok(norm(&g.map(D::re)) < tol) -} - -#[cfg(test)] -mod test { - use super::*; - use crate::gc_pcsaft::eos::parameter::test::*; - use approx::assert_relative_eq; - use feos_core::EosUnit; - use ndarray::arr1; - use num_dual::Dual64; - use quantity::si::{METER, MOL, PASCAL}; - - #[test] - fn test_assoc_propanol() { - let parameters = propanol(); - let contrib = Association { - parameters: Rc::new(parameters), - }; - let temperature = 300.0; - let volume = METER - .powi(3) - .to_reduced(EosUnit::reference_volume()) - .unwrap(); - let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); - let state = StateHD::new( - Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), - arr1(&[Dual64::from_re(moles)]), - ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); - assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); - } - - #[test] - fn test_cross_assoc_propanol() { - let parameters = propanol(); - let contrib = CrossAssociation { - parameters: Rc::new(parameters), - max_iter: 50, - tol: 1e-10, - }; - let temperature = 300.0; - let volume = METER - .powi(3) - .to_reduced(EosUnit::reference_volume()) - .unwrap(); - let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); - let state = StateHD::new( - Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), - arr1(&[Dual64::from_re(moles)]), - ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); - assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); - } - - #[test] - fn test_cross_assoc_ethanol_propanol() { - let parameters = ethanol_propanol(false); - let contrib = CrossAssociation { - parameters: Rc::new(parameters), - max_iter: 50, - tol: 1e-10, - }; - let temperature = 300.0; - let volume = METER - .powi(3) - .to_reduced(EosUnit::reference_volume()) - .unwrap(); - let moles = (arr1(&[1.5, 2.5]) * MOL) - .to_reduced(EosUnit::reference_moles()) - .unwrap(); - let state = StateHD::new( - Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), - moles.mapv(Dual64::from_re), - ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); - assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10); - } -} diff --git a/src/gc_pcsaft/eos/hard_sphere.rs b/src/gc_pcsaft/eos/hard_sphere.rs deleted file mode 100644 index c0455c25c..000000000 --- a/src/gc_pcsaft/eos/hard_sphere.rs +++ /dev/null @@ -1,129 +0,0 @@ -use super::GcPcSaftEosParameters; -use feos_core::{HelmholtzEnergyDual, StateHD}; -use ndarray::*; -use num_dual::DualNum; -use std::f64::consts::FRAC_PI_6; -use std::fmt; -use std::rc::Rc; - -impl GcPcSaftEosParameters { - pub fn hs_diameter>(&self, temperature: D) -> Array1 { - let ti = temperature.recip() * -3.0; - Array::from_shape_fn(self.sigma.len(), |i| { - -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] - }) - } -} - -#[derive(Clone)] -pub struct HardSphere { - pub parameters: Rc, -} - -impl> HelmholtzEnergyDual for HardSphere { - fn helmholtz_energy(&self, state: &StateHD) -> D { - let diameter = self.parameters.hs_diameter(state.temperature); - let zeta = self - .parameters - .zeta(&diameter, &state.partial_density, [0, 1, 2, 3]); - let frac_1mz3 = -(zeta[3] - 1.0).recip(); - let zeta_23 = self.parameters.zeta_23(&diameter, &state.molefracs); - state.volume * 6.0 / std::f64::consts::PI - * (zeta[1] * zeta[2] * frac_1mz3 * 3.0 - + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23 - + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p()) - } -} - -impl fmt::Display for HardSphere { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Hard Sphere (GC)") - } -} - -impl GcPcSaftEosParameters { - pub(crate) fn zeta, const N: usize>( - &self, - diameter: &Array1, - partial_density: &Array1, - k: [i32; N], - ) -> [D; N] { - let mut zeta = [D::zero(); N]; - for i in 0..self.m.len() { - for (z, &k) in zeta.iter_mut().zip(k.iter()) { - *z += partial_density[self.component_index[i]] - * diameter[i].powi(k) - * (FRAC_PI_6 * self.m[i]); - } - } - - zeta - } - - fn zeta_23>(&self, diameter: &Array1, molefracs: &Array1) -> D { - let mut zeta: [D; 2] = [D::zero(); 2]; - for i in 0..self.m.len() { - for (k, z) in zeta.iter_mut().enumerate() { - *z += molefracs[self.component_index[i]] - * diameter[i].powi((k + 2) as i32) - * (FRAC_PI_6 * self.m[i]); - } - } - - zeta[0] / zeta[1] - } -} - -#[cfg(test)] -mod test { - use super::*; - use crate::gc_pcsaft::eos::parameter::test::*; - use approx::assert_relative_eq; - use feos_core::EosUnit; - use num_dual::Dual64; - use quantity::si::{METER, MOL, PASCAL}; - - #[test] - fn test_hs_propane() { - let parameters = propane(); - let contrib = HardSphere { - parameters: Rc::new(parameters), - }; - let temperature = 300.0; - let volume = METER - .powi(3) - .to_reduced(EosUnit::reference_volume()) - .unwrap(); - let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); - let state = StateHD::new( - Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), - arr1(&[Dual64::from_re(moles)]), - ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); - assert_relative_eq!(pressure, 1.5285037907989527 * PASCAL, max_relative = 1e-10); - } - - #[test] - fn test_hs_propanol() { - let parameters = propanol(); - let contrib = HardSphere { - parameters: Rc::new(parameters), - }; - let temperature = 300.0; - let volume = METER - .powi(3) - .to_reduced(EosUnit::reference_volume()) - .unwrap(); - let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); - let state = StateHD::new( - Dual64::from_re(temperature), - Dual64::from_re(volume).derive(), - arr1(&[Dual64::from_re(moles)]), - ); - let pressure = - -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); - assert_relative_eq!(pressure, 2.3168212018200243 * PASCAL, max_relative = 1e-10); - } -} diff --git a/src/pcsaft/eos/association.rs b/src/pcsaft/eos/association.rs deleted file mode 100644 index 655c1c369..000000000 --- a/src/pcsaft/eos/association.rs +++ /dev/null @@ -1,295 +0,0 @@ -use super::PcSaftParameters; -use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; -use feos_saft::HardSphereProperties; -use ndarray::*; -use num_dual::linalg::{norm, LU}; -use num_dual::*; -use std::f64::consts::{FRAC_PI_3, PI}; -use std::fmt; -use std::rc::Rc; - -pub struct Association { - pub parameters: Rc, -} - -pub struct CrossAssociation { - pub parameters: Rc, - pub max_iter: usize, - pub tol: f64, -} - -fn association_strength>( - p: &PcSaftParameters, - temperature: D, - r: &Array1, - n2: D, - n3i: D, - xi: D, - ai: usize, - aj: usize, -) -> D { - let k = r[ai] * r[aj] / (r[ai] + r[aj]) * (n2 * n3i * 2.0); - n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) - * (0.5 * (p.sigma[ai] + p.sigma[aj])).powi(3) - * p.kappa_aibj[(ai, aj)] - * (temperature.recip() * p.epsilon_k_aibj[(ai, aj)]).exp_m1() -} - -impl> HelmholtzEnergyDual for Association { - fn helmholtz_energy(&self, state: &StateHD) -> D { - let p = &self.parameters; - let a = p.assoc_comp[0]; - - // temperature dependent segment radius - let r = p.hs_diameter(state.temperature) * 0.5; - - // auxiliary variables - let n2 = (&state.partial_density * &p.m * &r * &r).sum() * 4.0 * PI; - let n3 = (&state.partial_density * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3; - let n3i = (-n3 + 1.0).recip(); - - // association strength - let deltarho = association_strength(p, state.temperature, &r, n2, n3i, D::one(), a, a) - * state.partial_density[a]; - - let na = p.na[a]; - let nb = p.nb[a]; - if nb > 0.0 { - // no cross association, two association sites - let xa = assoc_site_frac_ab(deltarho, na, nb); - let xb = (xa - 1.0) * (na / nb) + 1.0; - - state.moles[a] * ((xa.ln() - xa * 0.5 + 0.5) * na + (xb.ln() - xb * 0.5 + 0.5) * nb) - } else { - // no cross association, one association site - let xa = assoc_site_frac_a(deltarho, na); - - state.moles[a] * (xa.ln() - xa * 0.5 + 0.5) * na - } - } -} - -impl fmt::Display for Association { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Association") - } -} - -pub(crate) fn assoc_site_frac_ab>(deltarho: D, na: f64, nb: f64) -> D { - if deltarho.re() > f64::EPSILON.sqrt() { - (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt() - - (deltarho * (nb - na) + 1.0)) - / (deltarho * na * 2.0) - } else { - D::one() + deltarho * nb * (deltarho * (nb + na) - 1.0) - } -} - -pub(crate) fn assoc_site_frac_a>(deltarho: D, na: f64) -> D { - if deltarho.re() > f64::EPSILON.sqrt() { - ((deltarho * na * 4.0 + 1.0).powi(2) - 1.0).sqrt() / (deltarho * na * 2.0) - } else { - D::one() + deltarho * na * (deltarho * na * 2.0 - 1.0) - } -} - -impl + ScalarOperand> HelmholtzEnergyDual for CrossAssociation { - fn helmholtz_energy(&self, state: &StateHD) -> D { - let p = &self.parameters; - - // temperature dependent segment radius - let r = p.hs_diameter(state.temperature) * 0.5; - - // auxiliary variables - let n2 = (&state.partial_density * &p.m * &r * &r).sum() * 4.0 * PI; - let n3 = (&state.partial_density * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3; - let n3i = (-n3 + 1.0).recip(); - - // Helmholtz energy - helmholtz_energy_density_cross_association( - p, - state.temperature, - &state.partial_density, - &r, - n2, - n3i, - D::one(), - self.max_iter, - self.tol, - None, - ) - .unwrap() - * state.volume - } -} - -impl fmt::Display for CrossAssociation { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Cross-association") - } -} - -pub fn helmholtz_energy_density_cross_association + ScalarOperand>( - p: &PcSaftParameters, - temperature: D, - density: &ArrayBase, - r: &Array1, - n2: D, - n3i: D, - xi: D, - max_iter: usize, - tol: f64, - x0: Option<&mut Array1>, -) -> Result -where - S: Data, -{ - // check if density is close to 0 - if density.sum().re() < f64::EPSILON { - if let Some(x0) = x0 { - x0.fill(1.0); - } - return Ok(D::zero()); - } - - // association strength - let delta = Array::from_shape_fn((p.nassoc, p.nassoc), |(i, j)| { - association_strength( - p, - temperature, - r, - n2, - n3i, - xi, - p.assoc_comp[i], - p.assoc_comp[j], - ) - }); - - // extract parameters of associating components - let na = Array::from_shape_fn(p.nassoc, |i| p.na[p.assoc_comp[i]]); - let nb = Array::from_shape_fn(p.nassoc, |i| p.nb[p.assoc_comp[i]]); - let rho = Array::from_shape_fn(p.nassoc, |i| density[p.assoc_comp[i]]); - - // cross-association according to Michelsen2006 - // initialize monomer fraction - let mut x = match &x0 { - Some(x0) => (*x0).clone(), - None => Array::from_elem(2 * p.nassoc, 0.2), - }; - - for k in 0..max_iter { - if newton_step_cross_association::( - &mut x, - p, - &delta.map(D::re), - &na, - &nb, - &rho.map(D::re), - tol, - )? { - break; - } - if k == max_iter - 1 { - return Err(EosError::NotConverged("Cross association".into())); - } - } - - // calculate derivatives - let mut x_dual = x.mapv(D::from); - for _ in 0..D::NDERIV { - newton_step_cross_association(&mut x_dual, p, &delta, &na, &nb, &rho, tol)?; - } - - // save monomer fraction - if let Some(x0) = x0 { - *x0 = x; - } - - // Helmholtz energy density - let xa = x_dual.slice(s![..p.nassoc]); - let xb = x_dual.slice(s![p.nassoc..]); - let f = |x: D| x.ln() - x * 0.5 + 0.5; - Ok((rho * (xa.mapv(f) * na + xb.mapv(f) * nb)).sum()) -} - -fn newton_step_cross_association + ScalarOperand>( - x: &mut Array1, - p: &PcSaftParameters, - delta: &Array2, - na: &Array1, - nb: &Array1, - rho: &Array1, - tol: f64, -) -> Result { - // gradient - let mut g: Array1 = Array::zeros(2 * p.nassoc); - // Hessian - let mut h: Array2 = Array::zeros((2 * p.nassoc, 2 * p.nassoc)); - - // slice arrays - let (xa, xb) = x.multi_slice_mut((s![..p.nassoc], s![p.nassoc..])); - let (mut ga, mut gb) = g.multi_slice_mut((s![..p.nassoc], s![p.nassoc..])); - let (mut haa, mut hab, mut hba, mut hbb) = h.multi_slice_mut(( - s![..p.nassoc, ..p.nassoc], - s![..p.nassoc, p.nassoc..], - s![p.nassoc.., ..p.nassoc], - s![p.nassoc.., p.nassoc..], - )); - - // calculate gradients and approximate Hessian - for i in 0..p.nassoc { - let d = &delta.index_axis(Axis(0), i) * rho; - - let dnx = (&xb * nb * &d).sum() + 1.0; - ga[i] = xa[i].recip() - dnx; - hab.index_axis_mut(Axis(0), i).assign(&(&d * &(-nb))); - haa[(i, i)] = -dnx / xa[i]; - - let dnx = (&xa * na * &d).sum() + 1.0; - gb[i] = xb[i].recip() - dnx; - hba.index_axis_mut(Axis(0), i).assign(&(&d * &(-na))); - hbb[(i, i)] = -dnx / xb[i]; - } - - // Newton step - x.assign(&(&*x - &LU::new(h)?.solve(&g))); - - // check convergence - Ok(norm(&g.map(D::re)) < tol) -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::pcsaft::parameters::utils::water_parameters; - use approx::assert_relative_eq; - - #[test] - fn helmholtz_energy() { - let assoc = Association { - parameters: Rc::new(water_parameters()), - }; - let t = 350.0; - let v = 41.248289328513216; - let n = 1.23; - let s = StateHD::new(t, v, arr1(&[n])); - let a_rust = assoc.helmholtz_energy(&s) / n; - assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); - } - - #[test] - fn helmholtz_energy_cross() { - let assoc = CrossAssociation { - parameters: Rc::new(water_parameters()), - max_iter: 50, - tol: 1e-10, - }; - let t = 350.0; - let v = 41.248289328513216; - let n = 1.23; - let s = StateHD::new(t, v, arr1(&[n])); - let a_rust = assoc.helmholtz_energy(&s) / n; - assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); - } -} diff --git a/src/pcsaft/eos/hard_sphere.rs b/src/pcsaft/eos/hard_sphere.rs deleted file mode 100644 index c52114029..000000000 --- a/src/pcsaft/eos/hard_sphere.rs +++ /dev/null @@ -1,112 +0,0 @@ -use super::PcSaftParameters; -use feos_core::{HelmholtzEnergyDual, StateHD}; -use ndarray::*; -use num_dual::DualNum; -use std::fmt; -use std::rc::Rc; - -impl PcSaftParameters { - pub fn hs_diameter>(&self, temperature: D) -> Array1 { - let ti = temperature.recip() * -3.0; - Array::from_shape_fn(self.sigma.len(), |i| { - -((ti * self.epsilon_k[i]).exp() * 0.12 - 1.0) * self.sigma[i] - }) - } -} - -pub struct HardSphere { - pub parameters: Rc, -} - -impl> HelmholtzEnergyDual for HardSphere { - fn helmholtz_energy(&self, state: &StateHD) -> D { - let d = self.parameters.hs_diameter(state.temperature); - let zeta = zeta(&self.parameters.m, &state.partial_density, &d); - let frac_1mz3 = -(zeta[3] - 1.0).recip(); - let zeta_23 = zeta_23(&self.parameters.m, &state.molefracs, &d); - state.volume * 6.0 / std::f64::consts::PI - * (zeta[1] * zeta[2] * frac_1mz3 * 3.0 - + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23 - + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p()) - } -} - -impl fmt::Display for HardSphere { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Hard Sphere") - } -} - -pub fn zeta>( - m: &Array1, - partial_density: &Array1, - diameter: &Array1, -) -> [D; 4] { - let mut zeta: [D; 4] = [D::zero(), D::zero(), D::zero(), D::zero()]; - for i in 0..m.len() { - for (k, z) in zeta.iter_mut().enumerate() { - *z += partial_density[i] - * diameter[i].powi(k as i32) - * (std::f64::consts::PI / 6.0 * m[i]); - } - } - zeta -} - -pub fn zeta_23>(m: &Array1, molefracs: &Array1, diameter: &Array1) -> D { - let mut zeta: [D; 2] = [D::zero(), D::zero()]; - for i in 0..m.len() { - for (k, z) in zeta.iter_mut().enumerate() { - *z += molefracs[i] * diameter[i].powi((k + 2) as i32) * m[i]; - } - } - zeta[0] / zeta[1] -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::pcsaft::parameters::utils::{ - butane_parameters, propane_butane_parameters, propane_parameters, - }; - use approx::assert_relative_eq; - use ndarray::arr1; - - #[test] - fn helmholtz_energy() { - let hs = HardSphere { - parameters: propane_parameters(), - }; - let t = 250.0; - let v = 1000.0; - let n = 1.0; - let s = StateHD::new(t, v, arr1(&[n])); - let a_rust = hs.helmholtz_energy(&s); - assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10); - } - - #[test] - fn mix() { - let c1 = HardSphere { - parameters: propane_parameters(), - }; - let c2 = HardSphere { - parameters: butane_parameters(), - }; - let c12 = HardSphere { - parameters: propane_butane_parameters(), - }; - let t = 250.0; - let v = 2.5e28; - let n = 1.0; - let s = StateHD::new(t, v, arr1(&[n])); - let a1 = c1.helmholtz_energy(&s); - let a2 = c2.helmholtz_energy(&s); - let s1m = StateHD::new(t, v, arr1(&[n, 0.0])); - let a1m = c12.helmholtz_energy(&s1m); - let s2m = StateHD::new(t, v, arr1(&[0.0, n])); - let a2m = c12.helmholtz_energy(&s2m); - assert_relative_eq!(a1, a1m, epsilon = 1e-14); - assert_relative_eq!(a2, a2m, epsilon = 1e-14); - } -} From a0e8fa49880381284003df841034f713c844aef0 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 4 Jul 2022 10:19:49 +0200 Subject: [PATCH 5/9] add missing files, update CI, and add tests from previously deleted files --- .github/workflows/release_saft.yml | 24 ++ .github/workflows/test.yml | 2 +- feos-saft/src/association/dft.rs | 166 ++++++++++ feos-saft/src/association/mod.rs | 479 ++++++++++++++++++++++++++++ feos-saft/src/association/python.rs | 46 +++ feos-saft/src/hard_sphere/dft.rs | 346 ++++++++++++++++++++ feos-saft/src/hard_sphere/mod.rs | 166 ++++++++++ src/gc_pcsaft/eos/mod.rs | 118 ++++++- src/pcsaft/eos/mod.rs | 63 +++- tests/gc_pcsaft/dft.rs | 53 --- 10 files changed, 1401 insertions(+), 62 deletions(-) create mode 100644 .github/workflows/release_saft.yml create mode 100644 feos-saft/src/association/dft.rs create mode 100644 feos-saft/src/association/mod.rs create mode 100644 feos-saft/src/association/python.rs create mode 100644 feos-saft/src/hard_sphere/dft.rs create mode 100644 feos-saft/src/hard_sphere/mod.rs diff --git a/.github/workflows/release_saft.yml b/.github/workflows/release_saft.yml new file mode 100644 index 000000000..26b9b4916 --- /dev/null +++ b/.github/workflows/release_saft.yml @@ -0,0 +1,24 @@ +name: Release feos-saft + +on: + push: + tags: ["feos-saft-v*"] + +jobs: + release-crates-io: + name: Release crates.io + runs-on: ubuntu-latest + defaults: + run: + working-directory: ./feos-saft + steps: + - uses: actions/checkout@v2 + - uses: actions-rs/toolchain@v1 + with: + profile: minimal + toolchain: stable + override: true + - uses: katyo/publish-crates@v1 + with: + registry-token: ${{ secrets.CRATES_IO_TOKEN }} + path: './feos-saft' diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 2ddd0138d..4170e477e 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -15,7 +15,7 @@ jobs: strategy: fail-fast: false matrix: - crate: [feos-core, feos-dft] + crate: [feos-core, feos-dft, feos-saft] steps: - uses: actions/checkout@v3 - name: Build diff --git a/feos-saft/src/association/dft.rs b/feos-saft/src/association/dft.rs new file mode 100644 index 000000000..5969dc32d --- /dev/null +++ b/feos-saft/src/association/dft.rs @@ -0,0 +1,166 @@ +use super::*; +use crate::HardSphereProperties; +use feos_core::EosResult; +use feos_dft::{ + FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, +}; +use num_dual::DualNum; +use std::f64::consts::PI; +use std::ops::MulAssign; + +pub const N0_CUTOFF: f64 = 1e-9; + +impl FunctionalContributionDual for Association

+where + N: DualNum + ScalarOperand, + P: HardSphereProperties, +{ + fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { + let p: &P = &self.parameters; + let r = p.hs_diameter(temperature) * 0.5; + let [_, _, _, c3] = p.geometry_coefficients(temperature); + WeightFunctionInfo::new(p.component_index().into_owned(), false) + .add( + WeightFunction::new_scaled(r.clone(), WeightFunctionShape::Delta), + false, + ) + .add( + WeightFunction { + prefactor: c3.clone(), + kernel_radius: r.clone(), + shape: WeightFunctionShape::DeltaVec, + }, + false, + ) + .add( + WeightFunction { + prefactor: c3, + kernel_radius: r, + shape: WeightFunctionShape::Theta, + }, + true, + ) + } + + fn calculate_helmholtz_energy_density( + &self, + temperature: N, + weighted_densities: ArrayView2, + ) -> EosResult> { + let p: &P = &self.parameters; + + // number of segments + let n = self.association_parameters.component_index.len(); + + // number of associating segments + let nassoc = self.association_parameters.assoc_comp.len(); + + // number of dimensions + let dim = (weighted_densities.shape()[0] - 1) / n - 1; + + // weighted densities + let n0i = weighted_densities.slice_axis(Axis(0), Slice::new(0, Some(n as isize), 1)); + let n2vi: Vec<_> = (0..dim) + .map(|i| { + weighted_densities.slice_axis( + Axis(0), + Slice::new((n * (i + 1)) as isize, Some((n * (i + 2)) as isize), 1), + ) + }) + .collect(); + let n3 = weighted_densities.index_axis(Axis(0), n * (dim + 1)); + + // calculate rho0 (only associating segments) + let [_, _, c2, _] = p.geometry_coefficients(temperature); + let diameter = p.hs_diameter(temperature); + let mut n2i = n0i.to_owned(); + for (i, mut n2i) in n2i.outer_iter_mut().enumerate() { + n2i.mul_assign(diameter[i].powi(2) * c2[i] * PI); + } + let mut rho0: Array2 = (n2vi + .iter() + .fold(Array::zeros(n0i.raw_dim()), |acc, n2vi| acc + n2vi * n2vi) + / -(&n2i * &n2i) + + 1.0) + * n0i; + rho0.iter_mut().zip(&n0i).for_each(|(rho0, &n0i)| { + if n0i.re() < N0_CUTOFF { + *rho0 = n0i; + } + }); + let rho0 = Array2::from_shape_fn((nassoc, n3.len()), |(i, j)| { + rho0[(self.association_parameters.assoc_comp[i], j)] + }); + + // calculate xi + let n2v: Vec<_> = n2vi.iter().map(|n2vi| n2vi.sum_axis(Axis(0))).collect(); + let n2 = n2i.sum_axis(Axis(0)); + let mut xi = n2v + .iter() + .fold(Array::zeros(n2.raw_dim()), |acc, n2v| acc + n2v) + / -(&n2 * &n2) + + 1.0; + xi.iter_mut() + .zip(&n0i.sum_axis(Axis(0))) + .for_each(|(xi, &n0i)| { + if n0i.re() < N0_CUTOFF { + *xi = N::one(); + } + }); + + // auxiliary variables + let n3i = n3.mapv(|n3| (-n3 + 1.0).recip()); + + // only one associating component + if nassoc == 1 { + // association strength + let k = &n2 * &n3i * diameter[self.association_parameters.assoc_comp[0]] * 0.5; + let deltarho = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) + * ((temperature.recip() * self.association_parameters.epsilon_k_aibj[(0, 0)]) + .exp_m1() + * self.association_parameters.sigma3_kappa_aibj[(0, 0)]) + * rho0.index_axis(Axis(0), 0); + + let na = self.association_parameters.na[0]; + let nb = self.association_parameters.nb[0]; + let f = |x: N| x.ln() - x * 0.5 + 0.5; + if nb > 0.0 { + // no cross association, two association sites + let xa = deltarho.mapv(|d| Self::assoc_site_frac_ab(d, na, nb)); + let xb = (&xa - 1.0) * (na / nb) + 1.0; + Ok((xa.mapv(f) * na + xb.mapv(f) * nb) * rho0.index_axis(Axis(0), 0)) + } else { + // no cross association, one association site + let xa = deltarho.mapv(|d| Self::assoc_site_frac_a(d, na)); + + Ok(xa.mapv(f) * na * rho0.index_axis(Axis(0), 0)) + } + } else { + let mut x: Array1 = Array::from_elem(2 * nassoc, 0.2); + Ok(rho0 + .view() + .into_shape([nassoc, rho0.len() / nassoc]) + .unwrap() + .axis_iter(Axis(1)) + .zip(n2.iter()) + .zip(n3i.iter()) + .zip(xi.iter()) + .map(|(((rho0, &n2), &n3i), &xi)| { + self.helmholtz_energy_density_cross_association( + temperature, + &rho0, + &diameter, + n2, + n3i, + xi, + self.max_iter, + self.tol, + Some(&mut x), + ) + }) + .collect::, _>>()? + .into_shape(n2.raw_dim()) + .unwrap()) + } + } +} diff --git a/feos-saft/src/association/mod.rs b/feos-saft/src/association/mod.rs new file mode 100644 index 000000000..aff05baee --- /dev/null +++ b/feos-saft/src/association/mod.rs @@ -0,0 +1,479 @@ +//! Generic implementation of the SAFT association contribution +//! that can be used across models. +use crate::HardSphereProperties; +use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; +use ndarray::*; +use num_dual::linalg::{norm, LU}; +use num_dual::*; +use serde::{Deserialize, Serialize}; +use std::fmt; +use std::rc::Rc; + +#[cfg(feature = "dft")] +pub(crate) mod dft; +#[cfg(feature = "python")] +pub(crate) mod python; + +/// Pure component association parameters +#[derive(Serialize, Deserialize, Clone, Default)] +pub struct AssociationRecord { + /// Association volume parameter + pub kappa_ab: f64, + /// Association energy parameter in units of Kelvin + pub epsilon_k_ab: f64, + /// \# of association sites of type A + pub na: f64, + /// \# of association sites of type B + pub nb: f64, +} + +impl AssociationRecord { + pub fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64) -> Self { + Self { + kappa_ab, + epsilon_k_ab, + na, + nb, + } + } +} + +impl fmt::Display for AssociationRecord { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "AssociationRecord(kappa_ab={}", self.kappa_ab)?; + write!(f, ", epsilon_k_ab={}", self.epsilon_k_ab)?; + write!(f, ", na={}", self.na)?; + write!(f, ", nb={})", self.nb) + } +} + +/// Parameter set required for the SAFT association Helmoltz energy +/// contribution and functional. +#[derive(Clone)] +pub struct AssociationParameters { + component_index: Array1, + pub assoc_comp: Array1, + pub kappa_ab: Array1, + pub epsilon_k_ab: Array1, + pub sigma3_kappa_aibj: Array2, + pub epsilon_k_aibj: Array2, + pub na: Array1, + pub nb: Array1, +} + +impl AssociationParameters { + pub fn new( + records: &[Option], + sigma: &Array1, + component_index: Option<&Array1>, + ) -> Self { + let mut assoc_comp = Vec::new(); + let mut sigma_assoc = Vec::new(); + let mut kappa_ab = Vec::new(); + let mut epsilon_k_ab = Vec::new(); + let mut na = Vec::new(); + let mut nb = Vec::new(); + + for (i, record) in records.iter().enumerate() { + if let Some(record) = record.as_ref() { + assoc_comp.push(i); + sigma_assoc.push(sigma[i]); + kappa_ab.push(record.kappa_ab); + epsilon_k_ab.push(record.epsilon_k_ab); + na.push(record.na); + nb.push(record.nb); + } + } + + let sigma3_kappa_aibj = Array2::from_shape_fn([kappa_ab.len(); 2], |(i, j)| { + (sigma_assoc[i] * sigma_assoc[j]).powf(1.5) * (kappa_ab[i] * kappa_ab[j]).sqrt() + }); + let epsilon_k_aibj = Array2::from_shape_fn([epsilon_k_ab.len(); 2], |(i, j)| { + 0.5 * (epsilon_k_ab[i] + epsilon_k_ab[j]) + }); + + Self { + component_index: component_index + .cloned() + .unwrap_or_else(|| Array1::from_shape_fn(assoc_comp.len(), |i| i)), + assoc_comp: Array1::from_vec(assoc_comp), + kappa_ab: Array1::from_vec(kappa_ab), + epsilon_k_ab: Array1::from_vec(epsilon_k_ab), + sigma3_kappa_aibj, + epsilon_k_aibj, + na: Array1::from_vec(na), + nb: Array1::from_vec(nb), + } + } +} + +/// Implementation of the SAFT association Helmholtz energy +/// contribution and functional. +pub struct Association

{ + parameters: Rc

, + association_parameters: AssociationParameters, + max_iter: usize, + tol: f64, + force_cross_association: bool, +} + +impl Association

{ + pub fn new( + parameters: &Rc

, + association_parameters: &AssociationParameters, + max_iter: usize, + tol: f64, + ) -> Self { + Self { + parameters: parameters.clone(), + association_parameters: association_parameters.clone(), + max_iter, + tol, + force_cross_association: false, + } + } + + pub fn new_cross_association( + parameters: &Rc

, + association_parameters: &AssociationParameters, + max_iter: usize, + tol: f64, + ) -> Self { + let mut res = Self::new(parameters, association_parameters, max_iter, tol); + res.force_cross_association = true; + res + } + + fn association_strength>( + &self, + temperature: D, + diameter: &Array1, + n2: D, + n3i: D, + xi: D, + ) -> Array2 { + // Calculate association strength + let ac = &self.association_parameters.assoc_comp; + Array2::from_shape_fn([ac.len(); 2], |(i, j)| { + let k = diameter[ac[i]] * diameter[ac[j]] / (diameter[ac[i]] + diameter[ac[j]]) + * (n2 * n3i); + n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) + * self.association_parameters.sigma3_kappa_aibj[(i, j)] + * (temperature.recip() * self.association_parameters.epsilon_k_aibj[(i, j)]) + .exp_m1() + }) + } +} + +impl + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual + for Association

+{ + fn helmholtz_energy(&self, state: &StateHD) -> D { + let p: &P = &self.parameters; + + // temperature dependent segment diameter + let diameter = p.hs_diameter(state.temperature); + + // auxiliary variables + let [zeta2, n3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); + let n2 = zeta2 * 6.0; + let n3i = (-n3 + 1.0).recip(); + + if self.association_parameters.assoc_comp.len() > 1 || self.force_cross_association { + // extract densities of associating segments + let rho_assoc = self + .association_parameters + .assoc_comp + .mapv(|a| state.partial_density[self.association_parameters.component_index[a]]); + + // Helmholtz energy + self.helmholtz_energy_density_cross_association( + state.temperature, + &rho_assoc, + &diameter, + n2, + n3i, + D::one(), + self.max_iter, + self.tol, + None, + ) + .unwrap() + * state.volume + } else { + // association strength + let c = self.association_parameters.component_index + [self.association_parameters.assoc_comp[0]]; + let deltarho = + self.association_strength(state.temperature, &diameter, n2, n3i, D::one())[(0, 0)] + * state.partial_density[c]; + + let na = self.association_parameters.na[0]; + let nb = self.association_parameters.nb[0]; + if nb > 0.0 { + // no cross association, two association sites + let xa = Self::assoc_site_frac_ab(deltarho, na, nb); + let xb = (xa - 1.0) * (na / nb) + 1.0; + + state.moles[c] * ((xa.ln() - xa * 0.5 + 0.5) * na + (xb.ln() - xb * 0.5 + 0.5) * nb) + } else { + // no cross association, one association site + let xa = Self::assoc_site_frac_a(deltarho, na); + + state.moles[c] * (xa.ln() - xa * 0.5 + 0.5) * na + } + } + } +} + +impl

fmt::Display for Association

{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Association") + } +} + +impl Association

{ + pub fn assoc_site_frac_ab>(deltarho: D, na: f64, nb: f64) -> D { + if deltarho.re() > f64::EPSILON.sqrt() { + (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt() + - (deltarho * (nb - na) + 1.0)) + / (deltarho * na * 2.0) + } else { + D::one() + deltarho * nb * (deltarho * (nb + na) - 1.0) + } + } + + pub fn assoc_site_frac_a>(deltarho: D, na: f64) -> D { + if deltarho.re() > f64::EPSILON.sqrt() { + ((deltarho * na * 4.0 + 1.0).powi(2) - 1.0).sqrt() / (deltarho * na * 2.0) + } else { + D::one() + deltarho * na * (deltarho * na * 2.0 - 1.0) + } + } + + #[allow(clippy::too_many_arguments)] + fn helmholtz_energy_density_cross_association< + S: Data, + D: DualNum + ScalarOperand, + >( + &self, + temperature: D, + density: &ArrayBase, + diameter: &Array1, + n2: D, + n3i: D, + xi: D, + max_iter: usize, + tol: f64, + x0: Option<&mut Array1>, + ) -> Result { + // check if density is close to 0 + if density.sum().re() < f64::EPSILON { + if let Some(x0) = x0 { + x0.fill(1.0); + } + return Ok(D::zero()); + } + + let assoc_comp = &self.association_parameters.assoc_comp; + let nassoc = assoc_comp.len(); + + // association strength + let delta = self.association_strength(temperature, diameter, n2, n3i, xi); + + // extract parameters of associating components + let na = &self.association_parameters.na; + let nb = &self.association_parameters.nb; + + // cross-association according to Michelsen2006 + // initialize monomer fraction + let mut x = match &x0 { + Some(x0) => (*x0).clone(), + None => Array::from_elem(2 * nassoc, 0.2), + }; + + for k in 0..max_iter { + if Self::newton_step_cross_association::<_, f64>( + nassoc, + &mut x, + &delta.map(D::re), + na, + nb, + &density.map(D::re), + tol, + )? { + break; + } + if k == max_iter - 1 { + return Err(EosError::NotConverged("Cross association".into())); + } + } + + // calculate derivatives + let mut x_dual = x.mapv(D::from); + for _ in 0..D::NDERIV { + Self::newton_step_cross_association(nassoc, &mut x_dual, &delta, na, nb, density, tol)?; + } + + // save monomer fraction + if let Some(x0) = x0 { + *x0 = x; + } + + // Helmholtz energy density + let xa = x_dual.slice(s![..nassoc]); + let xb = x_dual.slice(s![nassoc..]); + let f = |x: D| x.ln() - x * 0.5 + 0.5; + Ok((density * (xa.mapv(f) * na + xb.mapv(f) * nb)).sum()) + } + + fn newton_step_cross_association, D: DualNum + ScalarOperand>( + nassoc: usize, + x: &mut Array1, + delta: &Array2, + na: &Array1, + nb: &Array1, + rho: &ArrayBase, + tol: f64, + ) -> Result { + // gradient + let mut g: Array1 = Array::zeros(2 * nassoc); + // Hessian + let mut h: Array2 = Array::zeros((2 * nassoc, 2 * nassoc)); + + // slice arrays + let (xa, xb) = x.multi_slice_mut((s![..nassoc], s![nassoc..])); + let (mut ga, mut gb) = g.multi_slice_mut((s![..nassoc], s![nassoc..])); + let (mut haa, mut hab, mut hba, mut hbb) = h.multi_slice_mut(( + s![..nassoc, ..nassoc], + s![..nassoc, nassoc..], + s![nassoc.., ..nassoc], + s![nassoc.., nassoc..], + )); + + // calculate gradients and approximate Hessian + for i in 0..nassoc { + let d = &delta.index_axis(Axis(0), i) * rho; + + let dnx = (&xb * nb * &d).sum() + 1.0; + ga[i] = xa[i].recip() - dnx; + hab.index_axis_mut(Axis(0), i).assign(&(&d * &(-nb))); + haa[(i, i)] = -dnx / xa[i]; + + let dnx = (&xa * na * &d).sum() + 1.0; + gb[i] = xb[i].recip() - dnx; + hba.index_axis_mut(Axis(0), i).assign(&(&d * &(-na))); + hbb[(i, i)] = -dnx / xb[i]; + } + + // Newton step + x.assign(&(&*x - &LU::new(h)?.solve(&g))); + + // check convergence + Ok(norm(&g.map(D::re)) < tol) + } +} + +#[cfg(test)] +#[cfg(feature = "pcsaft")] +mod tests_pcsaft { + use super::*; + use crate::pcsaft::parameters::utils::water_parameters; + use approx::assert_relative_eq; + + #[test] + fn helmholtz_energy() { + let assoc = Association::new(&Rc::new(water_parameters()), 50, 1e-10); + let t = 350.0; + let v = 41.248289328513216; + let n = 1.23; + let s = StateHD::new(t, v, arr1(&[n])); + let a_rust = assoc.helmholtz_energy(&s) / n; + assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); + } + + #[test] + fn helmholtz_energy_cross() { + let assoc = Association::new_cross_association(&Rc::new(water_parameters()), 50, 1e-10); + let t = 350.0; + let v = 41.248289328513216; + let n = 1.23; + let s = StateHD::new(t, v, arr1(&[n])); + let a_rust = assoc.helmholtz_energy(&s) / n; + assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); + } +} + +#[cfg(test)] +#[cfg(feature = "gc_pcsaft")] +mod tests_gc_pcsaft { + use super::*; + use crate::gc_pcsaft::eos::parameter::test::*; + use approx::assert_relative_eq; + use feos_core::EosUnit; + use ndarray::arr1; + use num_dual::Dual64; + use quantity::si::{METER, MOL, PASCAL}; + + #[test] + fn test_assoc_propanol() { + let parameters = propanol(); + let contrib = Association::new(&Rc::new(parameters), 50, 1e-10); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + arr1(&[Dual64::from_re(moles)]), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); + } + + #[test] + fn test_cross_assoc_propanol() { + let parameters = propanol(); + let contrib = Association::new_cross_association(&Rc::new(parameters), 50, 1e-10); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + arr1(&[Dual64::from_re(moles)]), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); + } + + #[test] + fn test_cross_assoc_ethanol_propanol() { + let parameters = ethanol_propanol(false); + let contrib = Association::new(&Rc::new(parameters), 50, 1e-10); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (arr1(&[1.5, 2.5]) * MOL) + .to_reduced(EosUnit::reference_moles()) + .unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + moles.mapv(Dual64::from_re), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10); + } +} diff --git a/feos-saft/src/association/python.rs b/feos-saft/src/association/python.rs new file mode 100644 index 000000000..5a06f184c --- /dev/null +++ b/feos-saft/src/association/python.rs @@ -0,0 +1,46 @@ +use super::AssociationRecord; +use feos_core::impl_json_handling; +use feos_core::parameter::ParameterError; +use pyo3::prelude::*; + +/// Create a set of PC-Saft parameters from records. +#[pyclass(name = "AssociationRecord", unsendable)] +#[pyo3( + text_signature = "(m, sigma, epsilon_k, mu=None, q=None, kappa_ab=None, epsilon_k_ab=None, na=None, nb=None, viscosity=None, diffusion=None, thermal_conductivity=None)" +)] +#[derive(Clone)] +pub struct PyAssociationRecord(pub AssociationRecord); + +#[pymethods] +impl PyAssociationRecord { + #[new] + fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64) -> Self { + Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb)) + } + + #[getter] + fn get_kappa_ab(&self) -> f64 { + self.0.kappa_ab + } + + #[getter] + fn get_epsilon_k_ab(&self) -> f64 { + self.0.epsilon_k_ab + } + + #[getter] + fn get_na(&self) -> f64 { + self.0.na + } + + #[getter] + fn get_nb(&self) -> f64 { + self.0.nb + } + + fn __repr__(&self) -> PyResult { + Ok(self.0.to_string()) + } +} + +impl_json_handling!(PyAssociationRecord); diff --git a/feos-saft/src/hard_sphere/dft.rs b/feos-saft/src/hard_sphere/dft.rs new file mode 100644 index 000000000..6e346b810 --- /dev/null +++ b/feos-saft/src/hard_sphere/dft.rs @@ -0,0 +1,346 @@ +//! Helmholtz energy functionals from fundamental measure theory. +use feos_core::EosResult; +use feos_dft::adsorption::FluidParameters; +use feos_dft::solvation::PairPotential; +use feos_dft::{ + FunctionalContribution, FunctionalContributionDual, HelmholtzEnergyFunctional, MoleculeShape, + WeightFunction, WeightFunctionInfo, WeightFunctionShape, DFT, +}; +use ndarray::*; +use num_dual::DualNum; +use std::f64::consts::PI; +use std::fmt; +use std::rc::Rc; + +use crate::{HardSphereProperties, MonomerShape}; + +const PI36M1: f64 = 1.0 / (36.0 * PI); +const N3_CUTOFF: f64 = 1e-5; + +/// Different versions of fundamental measure theory +#[derive(Clone, Copy)] +#[cfg_attr(feature = "python", pyo3::pyclass)] +pub enum FMTVersion { + /// White Bear ([Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313)) or modified ([Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)) fundamental measure theory + WhiteBear, + /// Scalar fundamental measure theory by [Kierlik and Rosinberg, 1990](https://doi.org/10.1103/PhysRevA.42.3382) + KierlikRosinberg, + /// Anti-symmetric White Bear fundamental measure theory ([Rosenfeld et al., 1997](https://doi.org/10.1103/PhysRevE.55.4245)) and SI of ([Kessler et al., 2021](https://doi.org/10.1016/j.micromeso.2021.111263)) + AntiSymWhiteBear, +} + +/// The [FunctionalContribution] for the hard sphere functional. +pub struct FMTContribution

{ + pub properties: Rc

, + version: FMTVersion, +} + +impl

Clone for FMTContribution

{ + fn clone(&self) -> Self { + Self { + properties: self.properties.clone(), + version: self.version, + } + } +} + +impl

FMTContribution

{ + pub fn new(properties: &Rc

, version: FMTVersion) -> Self { + Self { + properties: properties.clone(), + version, + } + } +} + +impl> FunctionalContributionDual + for FMTContribution

+{ + fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { + let r = self.properties.hs_diameter(temperature) * 0.5; + let [c0, c1, c2, c3] = self.properties.geometry_coefficients(temperature); + match (self.version, r.len()) { + (FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear, 1) => { + WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) + .extend( + vec![ + WeightFunctionShape::Delta, + WeightFunctionShape::Theta, + WeightFunctionShape::DeltaVec, + ] + .into_iter() + .zip([c2, c3.clone(), c3]) + .map(|(s, c)| WeightFunction { + prefactor: c, + kernel_radius: r.clone(), + shape: s, + }) + .collect(), + false, + ) + } + (FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear, _) => { + WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) + .add( + WeightFunction { + prefactor: Zip::from(&c0) + .and(&r) + .map_collect(|&c, &r| r.powi(-2) * c / (4.0 * PI)), + kernel_radius: r.clone(), + shape: WeightFunctionShape::Delta, + }, + true, + ) + .add( + WeightFunction { + prefactor: Zip::from(&c1) + .and(&r) + .map_collect(|&c, &r| r.recip() * c / (4.0 * PI)), + kernel_radius: r.clone(), + shape: WeightFunctionShape::Delta, + }, + true, + ) + .add( + WeightFunction { + prefactor: c2, + kernel_radius: r.clone(), + shape: WeightFunctionShape::Delta, + }, + true, + ) + .add( + WeightFunction { + prefactor: c3.clone(), + kernel_radius: r.clone(), + shape: WeightFunctionShape::Theta, + }, + true, + ) + .add( + WeightFunction { + prefactor: Zip::from(&c3) + .and(&r) + .map_collect(|&c, &r| r.recip() * c / (4.0 * PI)), + kernel_radius: r.clone(), + shape: WeightFunctionShape::DeltaVec, + }, + true, + ) + .add( + WeightFunction { + prefactor: c3, + kernel_radius: r, + shape: WeightFunctionShape::DeltaVec, + }, + true, + ) + } + (FMTVersion::KierlikRosinberg, _) => { + WeightFunctionInfo::new(self.properties.component_index().into_owned(), false) + .extend( + vec![ + WeightFunctionShape::KR0, + WeightFunctionShape::KR1, + WeightFunctionShape::Delta, + WeightFunctionShape::Theta, + ] + .into_iter() + .zip(self.properties.geometry_coefficients(temperature)) + .map(|(s, c)| WeightFunction { + prefactor: c, + kernel_radius: r.clone(), + shape: s, + }) + .collect(), + true, + ) + } + } + } + + fn calculate_helmholtz_energy_density( + &self, + temperature: N, + weighted_densities: ArrayView2, + ) -> EosResult> { + let pure_component_weighted_densities = matches!( + self.version, + FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear + ) && self.properties.component_index().len() == 1; + + // scalar weighted densities + let (n2, n3) = if pure_component_weighted_densities { + ( + weighted_densities.index_axis(Axis(0), 0), + weighted_densities.index_axis(Axis(0), 1), + ) + } else { + ( + weighted_densities.index_axis(Axis(0), 2), + weighted_densities.index_axis(Axis(0), 3), + ) + }; + + let (n0, n1) = if pure_component_weighted_densities { + let r = self.properties.hs_diameter(temperature)[0] * 0.5; + ( + n2.mapv(|n2| n2 / (r * r * 4.0 * PI)), + n2.mapv(|n2| n2 / (r * 4.0 * PI)), + ) + } else { + ( + weighted_densities.index_axis(Axis(0), 0).to_owned(), + weighted_densities.index_axis(Axis(0), 1).to_owned(), + ) + }; + + // vector weighted densities + let (n1n2, n2n2) = match self.version { + FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear => { + let (n1v, n2v) = if pure_component_weighted_densities { + let r = self.properties.hs_diameter(temperature)[0] * 0.5; + let n2v = weighted_densities.slice_axis(Axis(0), Slice::new(2, None, 1)); + (n2v.mapv(|n2v| n2v / (r * 4.0 * PI)), n2v) + } else { + let dim = ((weighted_densities.shape()[0] - 4) / 2) as isize; + ( + weighted_densities + .slice_axis(Axis(0), Slice::new(4, Some(4 + dim), 1)) + .to_owned(), + weighted_densities + .slice_axis(Axis(0), Slice::new(4 + dim, Some(4 + 2 * dim), 1)), + ) + }; + match self.version { + FMTVersion::WhiteBear => ( + &n1 * &n2 - (&n1v * &n2v).sum_axis(Axis(0)), + &n2 * &n2 - (&n2v * &n2v).sum_axis(Axis(0)) * 3.0, + ), + FMTVersion::AntiSymWhiteBear => { + let mut xi2 = (&n2v * &n2v).sum_axis(Axis(0)) / n2.map(|n| n.powi(2)); + xi2.iter_mut().for_each(|x| { + if x.re() > 1.0 { + *x = N::one() + } + }); + ( + &n1 * &n2 - (&n1v * &n2v).sum_axis(Axis(0)), + &n2 * &n2 * xi2.mapv(|x| (-x + 1.0).powi(3)), + ) + } + _ => unreachable!(), + } + } + FMTVersion::KierlikRosinberg => (&n1 * &n2, &n2 * &n2), + }; + + // auxiliary variables + let ln31 = n3.mapv(|n3| (-n3).ln_1p()); + let n3rec = n3.mapv(|n3| n3.recip()); + let n3m1 = n3.mapv(|n3| -n3 + 1.0); + let n3m1rec = n3m1.mapv(|n3m1| n3m1.recip()); + + // use Taylor expansion for f3 at low densities to avoid numerical issues + let mut f3 = (&n3m1 * &n3m1 * &ln31 + n3) * &n3rec * n3rec * &n3m1rec * &n3m1rec; + f3.iter_mut().zip(n3).for_each(|(f3, &n3)| { + if n3.re() < N3_CUTOFF { + *f3 = (((n3 * 35.0 / 6.0 + 4.8) * n3 + 3.75) * n3 + 8.0 / 3.0) * n3 + 1.5; + } + }); + Ok(-(&n0 * &ln31) + n1n2 * &n3m1rec + n2n2 * n2 * PI36M1 * f3) + } +} + +impl fmt::Display for FMTContribution

{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let ver = match self.version { + FMTVersion::WhiteBear => "WB", + FMTVersion::KierlikRosinberg => "KR", + FMTVersion::AntiSymWhiteBear => "AntiSymWB", + }; + write!(f, "FMT functional ({})", ver) + } +} + +struct HardSphereParameters { + sigma: Array1, +} + +impl HardSphereProperties for HardSphereParameters { + fn monomer_shape(&self, _: N) -> MonomerShape { + MonomerShape::Spherical(self.sigma.len()) + } + + fn hs_diameter>(&self, _: N) -> Array1 { + self.sigma.mapv(N::from) + } +} + +/// [HelmholtzEnergyFunctional] for hard sphere systems. +pub struct FMTFunctional { + properties: Rc, + contributions: Vec>, + version: FMTVersion, +} + +impl FMTFunctional { + pub fn new(sigma: &Array1, version: FMTVersion) -> DFT { + let properties = Rc::new(HardSphereParameters { + sigma: sigma.clone(), + }); + let contributions: Vec> = + vec![Box::new(FMTContribution::new(&properties, version))]; + (Self { + properties, + contributions, + version, + }) + .into() + } +} + +impl HelmholtzEnergyFunctional for FMTFunctional { + fn contributions(&self) -> &[Box] { + &self.contributions + } + + fn subset(&self, component_list: &[usize]) -> DFT { + let sigma = component_list + .iter() + .map(|&c| self.properties.sigma[c]) + .collect(); + Self::new(&sigma, self.version) + } + + fn compute_max_density(&self, moles: &Array1) -> f64 { + moles.sum() / (moles * &self.properties.sigma).sum() * 1.2 + } + + fn molecule_shape(&self) -> MoleculeShape { + MoleculeShape::Spherical(self.properties.sigma.len()) + } +} + +impl PairPotential for FMTFunctional { + fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Array2 { + let s = &self.properties.sigma; + Array::from_shape_fn((s.len(), r.len()), |(j, k)| { + if r[k] > 0.5 * (s[i] + s[j]) { + 0.0 + } else { + f64::INFINITY + } + }) + } +} + +impl FluidParameters for FMTFunctional { + fn epsilon_k_ff(&self) -> Array1 { + Array::zeros(self.properties.sigma.len()) + } + + fn sigma_ff(&self) -> &Array1 { + &self.properties.sigma + } +} diff --git a/feos-saft/src/hard_sphere/mod.rs b/feos-saft/src/hard_sphere/mod.rs new file mode 100644 index 000000000..283abc729 --- /dev/null +++ b/feos-saft/src/hard_sphere/mod.rs @@ -0,0 +1,166 @@ +use feos_core::{HelmholtzEnergyDual, StateHD}; +use ndarray::*; +use num_dual::DualNum; +use std::borrow::Cow; +use std::f64::consts::FRAC_PI_6; +use std::fmt; +use std::rc::Rc; + +#[cfg(feature = "dft")] +mod dft; +#[cfg(feature = "dft")] +pub use dft::{FMTContribution, FMTFunctional, FMTVersion}; + +/// Different monomer shapes for FMT. +pub enum MonomerShape<'a, D> { + /// For spherical monomers, the number of components. + Spherical(usize), + /// For non-spherical molecules in a homosegmented approach, the + /// chain length parameter $m$. + NonSpherical(Array1), + /// For non-spherical molecules in a heterosegmented approach, + /// the geometry factors for every segment and the component + /// index for every segment. + Heterosegmented([Array1; 4], &'a Array1), +} + +/// Properties of (generalized) hard sphere systems. +pub trait HardSphereProperties { + fn monomer_shape>(&self, temperature: D) -> MonomerShape; + fn hs_diameter>(&self, temperature: D) -> Array1; + + fn component_index(&self) -> Cow> { + match self.monomer_shape(1.0) { + MonomerShape::Spherical(n) => Cow::Owned(Array1::from_shape_fn(n, |i| i)), + MonomerShape::NonSpherical(m) => Cow::Owned(Array1::from_shape_fn(m.len(), |i| i)), + MonomerShape::Heterosegmented(_, component_index) => Cow::Borrowed(component_index), + } + } + + fn geometry_coefficients>(&self, temperature: D) -> [Array1; 4] { + match self.monomer_shape(temperature) { + MonomerShape::Spherical(n) => { + let m = Array1::ones(n); + [m.clone(), m.clone(), m.clone(), m] + } + MonomerShape::NonSpherical(m) => [m.clone(), m.clone(), m.clone(), m], + MonomerShape::Heterosegmented(g, _) => g, + } + } + + fn zeta, const N: usize>( + &self, + temperature: D, + partial_density: &Array1, + k: [i32; N], + ) -> [D; N] { + let component_index = self.component_index(); + let geometry_coefficients = self.geometry_coefficients(temperature); + let diameter = self.hs_diameter(temperature); + let mut zeta = [D::zero(); N]; + for i in 0..diameter.len() { + for (z, &k) in zeta.iter_mut().zip(k.iter()) { + *z += partial_density[component_index[i]] + * diameter[i].powi(k) + * (geometry_coefficients[k as usize][i] * FRAC_PI_6); + } + } + + zeta + } + + fn zeta_23>(&self, temperature: D, molefracs: &Array1) -> D { + let component_index = self.component_index(); + let geometry_coefficients = self.geometry_coefficients(temperature); + let diameter = self.hs_diameter(temperature); + let mut zeta: [D; 2] = [D::zero(); 2]; + for i in 0..diameter.len() { + for (k, z) in zeta.iter_mut().enumerate() { + *z += molefracs[component_index[i]] + * diameter[i].powi((k + 2) as i32) + * (geometry_coefficients[k + 2][i] * FRAC_PI_6); + } + } + + zeta[0] / zeta[1] + } +} + +pub struct HardSphere

{ + parameters: Rc

, +} + +impl

HardSphere

{ + pub fn new(parameters: &Rc

) -> Self { + Self { + parameters: parameters.clone(), + } + } +} + +impl, P: HardSphereProperties> HelmholtzEnergyDual for HardSphere

{ + fn helmholtz_energy(&self, state: &StateHD) -> D { + let p = &self.parameters; + let zeta = p.zeta(state.temperature, &state.partial_density, [0, 1, 2, 3]); + let frac_1mz3 = -(zeta[3] - 1.0).recip(); + let zeta_23 = p.zeta_23(state.temperature, &state.molefracs); + state.volume * 6.0 / std::f64::consts::PI + * (zeta[1] * zeta[2] * frac_1mz3 * 3.0 + + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23 + + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p()) + } +} + +impl

fmt::Display for HardSphere

{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Hard Sphere") + } +} + +// #[cfg(test)] +// mod tests { +// use super::*; +// use crate::pcsaft::parameters::utils::{ +// butane_parameters, propane_butane_parameters, propane_parameters, +// }; +// use approx::assert_relative_eq; +// use ndarray::arr1; + +// #[test] +// fn helmholtz_energy() { +// let hs = HardSphere { +// parameters: propane_parameters(), +// }; +// let t = 250.0; +// let v = 1000.0; +// let n = 1.0; +// let s = StateHD::new(t, v, arr1(&[n])); +// let a_rust = hs.helmholtz_energy(&s); +// assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10); +// } + +// #[test] +// fn mix() { +// let c1 = HardSphere { +// parameters: propane_parameters(), +// }; +// let c2 = HardSphere { +// parameters: butane_parameters(), +// }; +// let c12 = HardSphere { +// parameters: propane_butane_parameters(), +// }; +// let t = 250.0; +// let v = 2.5e28; +// let n = 1.0; +// let s = StateHD::new(t, v, arr1(&[n])); +// let a1 = c1.helmholtz_energy(&s); +// let a2 = c2.helmholtz_energy(&s); +// let s1m = StateHD::new(t, v, arr1(&[n, 0.0])); +// let a1m = c12.helmholtz_energy(&s1m); +// let s2m = StateHD::new(t, v, arr1(&[0.0, n])); +// let a2m = c12.helmholtz_energy(&s2m); +// assert_relative_eq!(a1, a1m, epsilon = 1e-14); +// assert_relative_eq!(a2, a2m, epsilon = 1e-14); +// } +// } diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index edce4c102..05df696a2 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -52,9 +52,7 @@ impl GcPcSaft { pub fn with_options(parameters: Rc, options: GcPcSaftOptions) -> Self { let mut contributions: Vec> = Vec::with_capacity(7); - contributions.push(Box::new(HardSphere { - parameters: parameters.clone(), - })); + contributions.push(Box::new(HardSphere::new(¶meters))); contributions.push(Box::new(HardChain { parameters: parameters.clone(), })); @@ -117,3 +115,117 @@ impl MolarWeight for GcPcSaft { self.parameters.molarweight.clone() * GRAM / MOL } } + +#[cfg(test)] +mod test { + use super::*; + use crate::gc_pcsaft::eos::parameter::test::*; + use approx::assert_relative_eq; + use feos_core::{EosUnit, HelmholtzEnergyDual, StateHD}; + use ndarray::arr1; + use num_dual::Dual64; + use quantity::si::{METER, MOL, PASCAL}; + + #[test] + fn hs_propane() { + let parameters = propane(); + let contrib = HardSphere::new(&Rc::new(parameters)); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + arr1(&[Dual64::from_re(moles)]), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, 1.5285037907989527 * PASCAL, max_relative = 1e-10); + } + + #[test] + fn hs_propanol() { + let parameters = propanol(); + let contrib = HardSphere::new(&Rc::new(parameters)); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + arr1(&[Dual64::from_re(moles)]), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, 2.3168212018200243 * PASCAL, max_relative = 1e-10); + } + + #[test] + fn assoc_propanol() { + let parameters = Rc::new(propanol()); + let contrib = Association::new(¶meters, ¶meters.association, 50, 1e-10); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + arr1(&[Dual64::from_re(moles)]), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); + } + + #[test] + fn cross_assoc_propanol() { + let parameters = Rc::new(propanol()); + let contrib = + Association::new_cross_association(¶meters, ¶meters.association, 50, 1e-10); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (1.5 * MOL).to_reduced(EosUnit::reference_moles()).unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + arr1(&[Dual64::from_re(moles)]), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10); + } + + #[test] + fn cross_assoc_ethanol_propanol() { + let parameters = Rc::new(ethanol_propanol(false)); + let contrib = Association::new(¶meters, ¶meters.association, 50, 1e-10); + let temperature = 300.0; + let volume = METER + .powi(3) + .to_reduced(EosUnit::reference_volume()) + .unwrap(); + let moles = (arr1(&[1.5, 2.5]) * MOL) + .to_reduced(EosUnit::reference_moles()) + .unwrap(); + let state = StateHD::new( + Dual64::from_re(temperature), + Dual64::from_re(volume).derive(), + moles.mapv(Dual64::from_re), + ); + let pressure = + -contrib.helmholtz_energy(&state).eps[0] * temperature * EosUnit::reference_pressure(); + assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10); + } +} diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index d211c6500..a8b36d1a8 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -62,9 +62,7 @@ impl PcSaft { pub fn with_options(parameters: Rc, options: PcSaftOptions) -> Self { let mut contributions: Vec> = Vec::with_capacity(7); - contributions.push(Box::new(HardSphere { - parameters: parameters.clone(), - })); + contributions.push(Box::new(HardSphere::new(¶meters))); contributions.push(Box::new(HardChain { parameters: parameters.clone(), })); @@ -342,10 +340,10 @@ impl EntropyScaling for PcSaft { mod tests { use super::*; use crate::pcsaft::parameters::utils::{ - butane_parameters, propane_butane_parameters, propane_parameters, + butane_parameters, propane_butane_parameters, propane_parameters, water_parameters, }; use approx::assert_relative_eq; - use feos_core::{Contributions, DensityInitialization, PhaseEquilibrium, State}; + use feos_core::*; use ndarray::arr1; use quantity::si::{BAR, KELVIN, METER, PASCAL, RGAS, SECOND}; @@ -381,6 +379,61 @@ mod tests { ); } + #[test] + fn hard_sphere() { + let hs = HardSphere::new(&propane_parameters()); + let t = 250.0; + let v = 1000.0; + let n = 1.0; + let s = StateHD::new(t, v, arr1(&[n])); + let a_rust = hs.helmholtz_energy(&s); + assert_relative_eq!(a_rust, 0.410610492598808, epsilon = 1e-10); + } + + #[test] + fn hard_sphere_mix() { + let c1 = HardSphere::new(&propane_parameters()); + let c2 = HardSphere::new(&butane_parameters()); + let c12 = HardSphere::new(&propane_butane_parameters()); + let t = 250.0; + let v = 2.5e28; + let n = 1.0; + let s = StateHD::new(t, v, arr1(&[n])); + let a1 = c1.helmholtz_energy(&s); + let a2 = c2.helmholtz_energy(&s); + let s1m = StateHD::new(t, v, arr1(&[n, 0.0])); + let a1m = c12.helmholtz_energy(&s1m); + let s2m = StateHD::new(t, v, arr1(&[0.0, n])); + let a2m = c12.helmholtz_energy(&s2m); + assert_relative_eq!(a1, a1m, epsilon = 1e-14); + assert_relative_eq!(a2, a2m, epsilon = 1e-14); + } + + #[test] + fn association() { + let parameters = Rc::new(water_parameters()); + let assoc = Association::new(¶meters, ¶meters.association, 50, 1e-10); + let t = 350.0; + let v = 41.248289328513216; + let n = 1.23; + let s = StateHD::new(t, v, arr1(&[n])); + let a_rust = assoc.helmholtz_energy(&s) / n; + assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); + } + + #[test] + fn cross_association() { + let parameters = Rc::new(water_parameters()); + let assoc = + Association::new_cross_association(¶meters, ¶meters.association, 50, 1e-10); + let t = 350.0; + let v = 41.248289328513216; + let n = 1.23; + let s = StateHD::new(t, v, arr1(&[n])); + let a_rust = assoc.helmholtz_energy(&s) / n; + assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10); + } + #[test] fn new_tpn() { let e = Rc::new(PcSaft::new(propane_parameters())); diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index 69a865dd2..a09168eaa 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -169,59 +169,6 @@ fn test_bulk_association() -> Result<(), Box> { println!("{s:26}: {x:21.16}"); } assert_relative_eq!(p_eos[4].1, p_func[4].1, max_relative = 1e-14); - // println!("{:29}: {}", p_eos[0].0, p_eos[0].1); - // println!("{:29}: {}", p_func[0].0, p_func[0].1); - // println!(); - // println!("{:29}: {}", p_eos[1].0, p_eos[1].1); - // println!("{:29}: {}", p_func[1].0, p_func[1].1); - // println!(); - // println!("{:29}: {}", p_eos[2].0, p_eos[2].1); - // println!("{:29}: {}", p_func[2].0, p_func[2].1); - // println!(); - // println!("{:29}: {}", p_eos[3].0, p_eos[3].1); - // println!("{:29}: {}", p_func[3].0, p_func[3].1); - - // assert_relative_eq!( - // p_eos[0].1, - // 1.2471689792172869 * MEGA * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - // assert_relative_eq!( - // p_eos[1].1, - // 280.0635060891395938 * KILO * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - // assert_relative_eq!( - // p_eos[2].1, - // -141.9023918353318550 * KILO * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - // assert_relative_eq!( - // p_eos[3].1, - // -763.2289230004602132 * KILO * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - - // assert_relative_eq!( - // p_func[0].1, - // 1.2471689792172869 * MEGA * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - // assert_relative_eq!( - // p_func[1].1, - // 280.0635060891395938 * KILO * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - // assert_relative_eq!( - // p_func[2].1, - // -141.9023918353318550 * KILO * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); - // assert_relative_eq!( - // p_func[3].1, - // -763.2289230004602132 * KILO * PASCAL * KB / KB_old, - // max_relative = 1e-14, - // ); Ok(()) } From 644f3f3f4610fd293dd0035b0217213e51bd6168 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 4 Jul 2022 10:25:05 +0200 Subject: [PATCH 6/9] remove default features --- Cargo.toml | 2 +- feos-saft/Cargo.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index f0454b888..b0765b41b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,7 +47,7 @@ optional = true approx = "0.4" [features] -default = ["pcsaft", "dft", "gc_pcsaft"] +default = [] dft = ["feos-dft", "petgraph", "feos-saft?/dft"] estimator = [] pcsaft = ["feos-saft/association"] diff --git a/feos-saft/Cargo.toml b/feos-saft/Cargo.toml index fdc271d69..fdbf11f58 100644 --- a/feos-saft/Cargo.toml +++ b/feos-saft/Cargo.toml @@ -25,7 +25,7 @@ pyo3 = { version = "0.16", optional = true } serde = "1.0" [features] -default = ["association"] +default = [] association = [] dft = ["feos-dft"] python = ["pyo3"] \ No newline at end of file From ef6e28edc2b3df792cdb9d74ccca90255dd8feeb Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 4 Jul 2022 10:51:06 +0200 Subject: [PATCH 7/9] fix pets and python bindings --- Cargo.toml | 2 +- feos-saft/Cargo.toml | 1 + feos-saft/src/association/mod.rs | 6 ++++-- feos-saft/src/lib.rs | 2 ++ src/gc_pcsaft/python/mod.rs | 3 ++- src/pcsaft/python.rs | 3 ++- src/pets/dft/mod.rs | 12 +++--------- src/pets/dft/pure_pets_functional.rs | 2 +- tests/gc_pcsaft/dft.rs | 3 +-- 9 files changed, 17 insertions(+), 17 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index b0765b41b..d2d678801 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -53,6 +53,6 @@ estimator = [] pcsaft = ["feos-saft/association"] gc_pcsaft = ["feos-saft/association"] uvtheory = ["lazy_static"] -pets = [] +pets = ["feos-saft"] python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python", "feos-saft?/python"] all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets"] diff --git a/feos-saft/Cargo.toml b/feos-saft/Cargo.toml index fdbf11f58..cbb54f9a2 100644 --- a/feos-saft/Cargo.toml +++ b/feos-saft/Cargo.toml @@ -23,6 +23,7 @@ feos-dft = { version = "0.2", path = "../feos-dft", optional = true } ndarray = { version = "0.15", features = ["serde"] } pyo3 = { version = "0.16", optional = true } serde = "1.0" +serde_json = "1.0" [features] default = [] diff --git a/feos-saft/src/association/mod.rs b/feos-saft/src/association/mod.rs index aff05baee..aabfeffb0 100644 --- a/feos-saft/src/association/mod.rs +++ b/feos-saft/src/association/mod.rs @@ -10,9 +10,11 @@ use std::fmt; use std::rc::Rc; #[cfg(feature = "dft")] -pub(crate) mod dft; +mod dft; #[cfg(feature = "python")] -pub(crate) mod python; +mod python; +#[cfg(feature = "python")] +pub use python::PyAssociationRecord; /// Pure component association parameters #[derive(Serialize, Deserialize, Clone, Default)] diff --git a/feos-saft/src/lib.rs b/feos-saft/src/lib.rs index fe2d0d59f..c8cdd98fd 100644 --- a/feos-saft/src/lib.rs +++ b/feos-saft/src/lib.rs @@ -5,5 +5,7 @@ pub use hard_sphere::{HardSphere, HardSphereProperties, MonomerShape}; #[cfg(feature = "association")] mod association; +#[cfg(all(feature = "association", feature = "python"))] +pub use association::PyAssociationRecord; #[cfg(feature = "association")] pub use association::{Association, AssociationParameters, AssociationRecord}; diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs index 3b874a6a2..5164eed57 100644 --- a/src/gc_pcsaft/python/mod.rs +++ b/src/gc_pcsaft/python/mod.rs @@ -2,7 +2,6 @@ use super::dft::GcPcSaftFunctionalParameters; use super::eos::GcPcSaftEosParameters; use super::record::GcPcSaftRecord; -use crate::association::python::PyAssociationRecord; use feos_core::joback::JobackRecord; use feos_core::parameter::{ BinaryRecord, IdentifierOption, ParameterError, ParameterHetero, SegmentRecord, @@ -10,6 +9,7 @@ use feos_core::parameter::{ use feos_core::python::joback::PyJobackRecord; use feos_core::python::parameter::{PyBinarySegmentRecord, PyChemicalRecord, PyIdentifier}; use feos_core::{impl_json_handling, impl_parameter_from_segments, impl_segment_record}; +use feos_saft::PyAssociationRecord; #[cfg(feature = "dft")] use numpy::{PyArray2, ToPyArray}; use pyo3::prelude::*; @@ -151,6 +151,7 @@ pub fn gc_pcsaft(_py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index 6f70a484d..66d83a46e 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -1,5 +1,4 @@ use super::parameters::{PcSaftBinaryRecord, PcSaftParameters, PcSaftRecord}; -use crate::association::python::PyAssociationRecord; use feos_core::joback::JobackRecord; use feos_core::parameter::{ BinaryRecord, Identifier, IdentifierOption, Parameter, ParameterError, PureRecord, @@ -8,6 +7,7 @@ use feos_core::parameter::{ use feos_core::python::joback::PyJobackRecord; use feos_core::python::parameter::*; use feos_core::*; +use feos_saft::PyAssociationRecord; use numpy::{PyArray2, PyReadonlyArray2, ToPyArray}; use pyo3::exceptions::PyTypeError; use pyo3::prelude::*; @@ -161,6 +161,7 @@ pub fn pcsaft(_py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/src/pets/dft/mod.rs b/src/pets/dft/mod.rs index 222c98909..f70ea62a9 100644 --- a/src/pets/dft/mod.rs +++ b/src/pets/dft/mod.rs @@ -5,12 +5,10 @@ use feos_core::joback::Joback; use feos_core::parameter::Parameter; use feos_core::{IdealGasContribution, MolarWeight}; use feos_dft::adsorption::FluidParameters; -use feos_dft::fundamental_measure_theory::{ - FMTContribution, FMTProperties, FMTVersion, MonomerShape, -}; use feos_dft::solvation::PairPotential; use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use ndarray::{Array, Array1, Array2}; +use feos_saft::{FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; +use ndarray::{Array1, Array2}; use num_dual::DualNum; use pure_pets_functional::*; use quantity::si::*; @@ -117,11 +115,7 @@ impl MolarWeight for PetsFunctional { } } -impl FMTProperties for PetsParameters { - fn component_index(&self) -> Array1 { - Array::from_shape_fn(self.sigma.len(), |i| i) - } - +impl HardSphereProperties for PetsParameters { fn monomer_shape>(&self, _: N) -> MonomerShape { MonomerShape::Spherical(self.sigma.len()) } diff --git a/src/pets/dft/pure_pets_functional.rs b/src/pets/dft/pure_pets_functional.rs index 00b566112..c88fc9e9e 100644 --- a/src/pets/dft/pure_pets_functional.rs +++ b/src/pets/dft/pure_pets_functional.rs @@ -1,10 +1,10 @@ use crate::pets::eos::dispersion::{A, B}; use crate::pets::parameters::PetsParameters; use feos_core::{EosError, EosResult}; -use feos_dft::fundamental_measure_theory::FMTVersion; use feos_dft::{ FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, }; +use feos_saft::FMTVersion; use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index a09168eaa..cf0a77a98 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -3,10 +3,9 @@ use approx::assert_relative_eq; use feos::gc_pcsaft::{ GcPcSaft, GcPcSaftEosParameters, GcPcSaftFunctional, GcPcSaftFunctionalParameters, - GcPcSaftRecord, }; use feos_core::parameter::{ - self, ChemicalRecord, Identifier, IdentifierOption, ParameterHetero, SegmentRecord, + ChemicalRecord, Identifier, IdentifierOption, ParameterHetero, SegmentRecord, }; use feos_core::{PhaseEquilibrium, State, StateBuilder, Verbosity}; use feos_dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification}; From d7a217be67efd5cf45eff43e035f7a4e9c9778ed Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 4 Jul 2022 20:28:35 +0200 Subject: [PATCH 8/9] small but very important fix --- feos-saft/src/association/dft.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/feos-saft/src/association/dft.rs b/feos-saft/src/association/dft.rs index 5969dc32d..e143ac00b 100644 --- a/feos-saft/src/association/dft.rs +++ b/feos-saft/src/association/dft.rs @@ -16,7 +16,7 @@ where P: HardSphereProperties, { fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { - let p: &P = &self.parameters; + let p = &self.parameters; let r = p.hs_diameter(temperature) * 0.5; let [_, _, _, c3] = p.geometry_coefficients(temperature); WeightFunctionInfo::new(p.component_index().into_owned(), false) @@ -47,7 +47,7 @@ where temperature: N, weighted_densities: ArrayView2, ) -> EosResult> { - let p: &P = &self.parameters; + let p = &self.parameters; // number of segments let n = self.association_parameters.component_index.len(); @@ -97,7 +97,7 @@ where let n2 = n2i.sum_axis(Axis(0)); let mut xi = n2v .iter() - .fold(Array::zeros(n2.raw_dim()), |acc, n2v| acc + n2v) + .fold(Array::zeros(n2.raw_dim()), |acc, n2v| acc + n2v * n2v) / -(&n2 * &n2) + 1.0; xi.iter_mut() From 2ac40f13914127653fa7a9e25baebf8a5afb1a28 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Wed, 13 Jul 2022 09:14:56 +0200 Subject: [PATCH 9/9] update changelog and fix bug in wheels --- feos-dft/CHANGELOG.md | 3 +++ src/pcsaft/dft/mod.rs | 2 +- src/pets/dft/mod.rs | 2 +- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index c3537a1ab..95a622910 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -13,6 +13,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Made FMT functional more flexible w.r.t. the shape of the weight functions. [#31](https://github.com/feos-org/feos-dft/pull/31) - Changed interface of `PairCorrelationFunction` to facilitate the calculation of pair correlation functions in mixtures. [#29](https://github.com/feos-org/feos-dft/pull/29) +### Removed +- Moved the implementation of fundamental measure theory to the new `feos-saft` crate. [#27](https://github.com/feos-org/feos/pull/27) + ## [0.2.0] - 2022-04-12 ### Added - Added `grand_potential_density` getter for DFT profiles in Python. [#22](https://github.com/feos-org/feos-dft/pull/22) diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index 257597dfd..2a520ee8e 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -150,7 +150,7 @@ impl PairPotential for PcSaftFunctional { fn pair_potential(&self, i: usize, r: &Array1, _: f64) -> Array2 { let sigma_ij = &self.parameters.sigma_ij; let eps_ij_4 = 4.0 * &self.parameters.epsilon_k_ij; - Array::from_shape_fn((self.parameters.m.len(), r.len()), |(j, k)| { + Array2::from_shape_fn((self.parameters.m.len(), r.len()), |(j, k)| { let att = (sigma_ij[[i, j]] / r[k]).powi(6); eps_ij_4[[i, j]] * att * (att - 1.0) }) diff --git a/src/pets/dft/mod.rs b/src/pets/dft/mod.rs index f70ea62a9..9ed66588e 100644 --- a/src/pets/dft/mod.rs +++ b/src/pets/dft/mod.rs @@ -140,7 +140,7 @@ impl PairPotential for PetsFunctional { let eps_ij_4 = 4.0 * self.parameters.epsilon_k_ij.clone(); let shift_ij = &eps_ij_4 * (2.5.powi(-12) - 2.5.powi(-6)); let rc_ij = 2.5 * &self.parameters.sigma_ij; - Array::from_shape_fn((self.parameters.sigma.len(), r.len()), |(j, k)| { + Array2::from_shape_fn((self.parameters.sigma.len(), r.len()), |(j, k)| { if r[k] > rc_ij[[i, j]] { 0.0 } else {