From fe17f689f8e757e08acc5e93b9f37daea36e1069 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 13 Mar 2025 11:47:09 +0100 Subject: [PATCH 1/3] Remove the `DFT` wrapper struct --- feos-core/src/python/phase_equilibria.rs | 2 +- feos-derive/src/dft.rs | 172 ++++---------- feos-derive/src/eos.rs | 210 ------------------ feos-derive/src/functional_contribution.rs | 6 +- feos-derive/src/lib.rs | 12 +- feos-derive/src/residual.rs | 60 +++-- feos-dft/src/adsorption/mod.rs | 18 +- feos-dft/src/adsorption/pore.rs | 33 ++- feos-dft/src/adsorption/pore2d.rs | 4 +- feos-dft/src/adsorption/pore3d.rs | 4 +- feos-dft/src/functional.rs | 112 ++-------- feos-dft/src/interface/mod.rs | 15 +- .../src/interface/surface_tension_diagram.rs | 8 +- feos-dft/src/lib.rs | 3 +- feos-dft/src/pdgt.rs | 11 +- feos-dft/src/profile/mod.rs | 10 +- feos-dft/src/solvation/pair_correlation.rs | 4 +- feos-dft/src/solvation/solvation_profile.rs | 4 +- src/eos.rs | 80 +++++-- src/functional.rs | 55 ----- src/gc_pcsaft/dft/mod.rs | 40 ++-- src/gc_pcsaft/micelles.rs | 6 +- src/hard_sphere/dft.rs | 29 ++- src/lib.rs | 4 - src/pcsaft/dft/mod.rs | 34 +-- src/pets/dft/mod.rs | 38 ++-- src/python/dft.rs | 125 +++++------ src/python/eos.rs | 2 +- src/python/mod.rs | 3 - src/saftvrqmie/dft/mod.rs | 34 +-- tests/pcsaft/dft.rs | 9 +- 31 files changed, 405 insertions(+), 742 deletions(-) delete mode 100644 feos-derive/src/eos.rs delete mode 100644 src/functional.rs diff --git a/feos-core/src/python/phase_equilibria.rs b/feos-core/src/python/phase_equilibria.rs index e9b5e00f2..fb761c2da 100644 --- a/feos-core/src/python/phase_equilibria.rs +++ b/feos-core/src/python/phase_equilibria.rs @@ -4,7 +4,7 @@ macro_rules! impl_phase_equilibrium { /// A thermodynamic two phase equilibrium state. #[pyclass(name = "PhaseEquilibrium")] #[derive(Clone)] - pub struct PyPhaseEquilibrium(PhaseEquilibrium<$eos, 2>); + pub struct PyPhaseEquilibrium(pub PhaseEquilibrium<$eos, 2>); #[pymethods] impl PyPhaseEquilibrium { diff --git a/feos-derive/src/dft.rs b/feos-derive/src/dft.rs index d71ec58ba..a63974b28 100644 --- a/feos-derive/src/dft.rs +++ b/feos-derive/src/dft.rs @@ -1,114 +1,45 @@ +use crate::{implement, OPT_IMPLS}; use quote::quote; use syn::DeriveInput; -use crate::implement; - -const OPT_IMPLS: [&str; 4] = [ - "bond_lengths", - "molar_weight", - "fluid_parameters", - "pair_potential", -]; - pub(crate) fn expand_helmholtz_energy_functional( input: DeriveInput, ) -> syn::Result { - let variants = match input.data { - syn::Data::Enum(syn::DataEnum { ref variants, .. }) => variants, - _ => panic!("this derive macro only works on enums"), + let syn::Data::Enum(syn::DataEnum { ref variants, .. }) = input.data else { + panic!("this derive macro only works on enums") }; - let from = impl_from(variants)?; - let functional = impl_helmholtz_energy_functional(variants)?; - let fluid_parameters = impl_fluid_parameters(variants)?; - let pair_potential = impl_pair_potential(variants)?; + let functional = impl_helmholtz_energy_functional(&input.ident, variants)?; + let fluid_parameters = impl_fluid_parameters(&input.ident, variants)?; + let pair_potential = impl_pair_potential(&input.ident, variants)?; Ok(quote! { - #from #functional #fluid_parameters #pair_potential }) } -// extract the variant name and the name of the functional, -// i.e. PcSaft(PcSaftFunctional) will return (PcSaft, PcSaftFunctional) -fn extract_names(variant: &syn::Variant) -> syn::Result<(&syn::Ident, &syn::Ident)> { - let name = &variant.ident; - let field = if let syn::Fields::Unnamed(syn::FieldsUnnamed { ref unnamed, .. }) = variant.fields - { - if unnamed.len() != 1 { - return Err(syn::Error::new_spanned( - unnamed, - "expected tuple struct with single HelmholtzFunctional as variant", - )); - } - &unnamed[0] - } else { - return Err(syn::Error::new_spanned( - name, - "expected variant with a HelmholtzFunctional as data", - )); - }; - - let inner = if let syn::Type::Path(syn::TypePath { ref path, .. }) = &field.ty { - path.get_ident() - } else { - None - } - .ok_or_else(|| syn::Error::new_spanned(field, "expected HelmholtzFunctional"))?; - Ok((name, inner)) -} - -fn impl_from( - variants: &syn::punctuated::Punctuated, -) -> syn::Result { - variants - .iter() - .map(|v| { - let (variant_name, functional_name) = extract_names(v)?; - Ok(quote! { - impl From<#functional_name> for FunctionalVariant { - fn from(f: #functional_name) -> Self { - Self::#variant_name(f) - } - } - }) - }) - .collect() -} - -fn impl_helmholtz_energy_functional( +pub(crate) fn impl_helmholtz_energy_functional( + ident: &syn::Ident, variants: &syn::punctuated::Punctuated, ) -> syn::Result { - let molecule_shape = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(functional) => functional.molecule_shape() - } - }); - let compute_max_density = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(functional) => functional.compute_max_density(moles) - } - }); - let contributions = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(functional) => Box::new(functional.contributions().map(FunctionalContributionVariant::from)) - } - }); - - let mut molar_weight = Vec::new(); - let mut has_molar_weight = Vec::new(); + let mut molecule_shape = Vec::new(); + let mut contributions = Vec::new(); for v in variants.iter() { - if implement("molar_weight", v, &OPT_IMPLS)? { - let name = &v.ident; - molar_weight.push(quote! { - Self::#name(functional) => functional.molar_weight() + let name = &v.ident; + if implement("functional", v, &OPT_IMPLS)? { + molecule_shape.push(quote! { + Self::#name(functional) => functional.molecule_shape() + }); + contributions.push(quote! { + Self::#name(functional) => Box::new(functional.contributions().map(FunctionalContributionVariant::from)) + }); + } else { + molecule_shape.push(quote! { + Self::#name(functional) => panic!("{} is not a Helmholtz energy functional!", stringify!(#name)) }); - has_molar_weight.push(quote! { - Self::#name(functional) => true + contributions.push(quote! { + Self::#name(functional) => panic!("{} is not a Helmholtz energy functional!", stringify!(#name)) }); } } @@ -124,45 +55,22 @@ fn impl_helmholtz_energy_functional( } Ok(quote! { - impl HelmholtzEnergyFunctional for FunctionalVariant { + impl HelmholtzEnergyFunctional for #ident { type Contribution = FunctionalContributionVariant; - fn molecule_shape(&self) -> MoleculeShape { + fn molecule_shape(&self) -> feos_dft::MoleculeShape { match self { #(#molecule_shape,)* } } - fn compute_max_density(&self, moles: &Array1) -> f64 { - match self { - #(#compute_max_density,)* - } - } fn contributions(&self) -> Box> { match self { #(#contributions,)* } } - fn bond_lengths + Copy>(&self, temperature: N) -> UnGraph<(), N> { + fn bond_lengths + Copy>(&self, temperature: N) -> petgraph::graph::UnGraph<(), N> { match self { #(#bond_lengths,)* - _ => Graph::with_capacity(0, 0), - } - } - } - - impl Molarweight for FunctionalVariant { - fn molar_weight(&self) -> MolarWeight> { - match self { - #(#molar_weight,)* - _ => unimplemented!() - } - } - } - - impl FunctionalVariant { - pub fn has_molar_weight(&self) -> bool { - match self { - #(#has_molar_weight,)* - _ => false, + _ => petgraph::Graph::with_capacity(0, 0), } } } @@ -170,35 +78,41 @@ fn impl_helmholtz_energy_functional( } fn impl_fluid_parameters( + ident: &syn::Ident, variants: &syn::punctuated::Punctuated, ) -> syn::Result { let mut epsilon_k_ff = Vec::new(); let mut sigma_ff = Vec::new(); for v in variants.iter() { + let name = &v.ident; if implement("fluid_parameters", v, &OPT_IMPLS)? { - let name = &v.ident; epsilon_k_ff.push(quote! { Self::#name(functional) => functional.epsilon_k_ff() }); sigma_ff.push(quote! { Self::#name(functional) => functional.sigma_ff() }); + } else { + epsilon_k_ff.push(quote! { + Self::#name(functional) => panic!("{} does not support the automatic calculation of external potentials!", stringify!(#name)) + }); + sigma_ff.push(quote! { + Self::#name(functional) => panic!("{} does not support the automatic calculation of external potentials!", stringify!(#name)) + }); } } Ok(quote! { - impl FluidParameters for FunctionalVariant { + impl feos_dft::adsorption::FluidParameters for #ident { fn epsilon_k_ff(&self) -> Array1 { match self { #(#epsilon_k_ff,)* - _ => unimplemented!() } } fn sigma_ff(&self) -> &Array1 { match self { #(#sigma_ff,)* - _ => unimplemented!() } } } @@ -206,24 +120,28 @@ fn impl_fluid_parameters( } fn impl_pair_potential( + ident: &syn::Ident, variants: &syn::punctuated::Punctuated, ) -> syn::Result { let mut pair_potential = Vec::new(); for v in variants.iter() { + let name = &v.ident; if implement("pair_potential", v, &OPT_IMPLS)? { - let name = &v.ident; pair_potential.push(quote! { Self::#name(functional) => functional.pair_potential(i, r, temperature) }); + } else { + pair_potential.push(quote! { + Self::#name(functional) => panic!("{} does not provide pair potentials!", stringify!(#name)) + }); } } Ok(quote! { - impl PairPotential for FunctionalVariant { - fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> Array2 { + impl feos_dft::solvation::PairPotential for #ident { + fn pair_potential(&self, i: usize, r: &Array1, temperature: f64) -> ndarray::Array2 { match self { #(#pair_potential,)* - _ => unimplemented!() } } } diff --git a/feos-derive/src/eos.rs b/feos-derive/src/eos.rs deleted file mode 100644 index 6dc4afea1..000000000 --- a/feos-derive/src/eos.rs +++ /dev/null @@ -1,210 +0,0 @@ -use super::implement; -use quote::quote; -use syn::DeriveInput; - -// possible additional traits to implement -const OPT_IMPLS: [&str; 2] = ["molar_weight", "entropy_scaling"]; - -pub(crate) fn expand_equation_of_state( - input: DeriveInput, -) -> syn::Result { - let variants = match input.data { - syn::Data::Enum(syn::DataEnum { ref variants, .. }) => variants, - _ => panic!("this derive macro only works on enums"), - }; - - let eos = impl_equation_of_state(variants); - let molar_weight = impl_molar_weight(variants)?; - let entropy_scaling = impl_entropy_scaling(variants)?; - Ok(quote! { - #eos - #molar_weight - #entropy_scaling - }) -} - -fn impl_equation_of_state( - variants: &syn::punctuated::Punctuated, -) -> proc_macro2::TokenStream { - let components = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(eos) => eos.components() - } - }); - let compute_max_density = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(eos) => eos.compute_max_density(moles) - } - }); - let subset = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(eos) => Self::#name(eos.subset(component_list)) - } - }); - let residual = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(eos) => eos.residual() - } - }); - let ideal_gas = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(eos) => eos.ideal_gas() - } - }); - - quote! { - impl EquationOfState for EosVariant { - fn components(&self) -> usize { - match self { - #(#components,)* - } - } - fn compute_max_density(&self, moles: &Array1) -> f64 { - match self { - #(#compute_max_density,)* - } - } - fn subset(&self, component_list: &[usize]) -> Self { - match self { - #(#subset,)* - } - } - fn residual(&self) -> &[Box] { - match self { - #(#residual,)* - } - } - fn ideal_gas(&self) -> &dyn IdealGasContribution { - match self { - #(#ideal_gas,)* - } - } - } - } -} - -fn impl_molar_weight( - variants: &syn::punctuated::Punctuated, -) -> syn::Result { - let mut molar_weight = Vec::new(); - - for v in variants.iter() { - if implement("molar_weight", v, &OPT_IMPLS)? { - let name = &v.ident; - molar_weight.push(quote! { - Self::#name(eos) => eos.molar_weight() - }); - } - } - Ok(quote! { - impl MolarWeight for EosVariant { - fn molar_weight(&self) -> SIArray1 { - match self { - #(#molar_weight,)* - _ => unimplemented!() - } - } - } - }) -} - -fn impl_entropy_scaling( - variants: &syn::punctuated::Punctuated, -) -> syn::Result { - let mut etar = Vec::new(); - let mut etac = Vec::new(); - let mut dr = Vec::new(); - let mut dc = Vec::new(); - let mut thcr = Vec::new(); - let mut thcc = Vec::new(); - - for v in variants.iter() { - if implement("entropy_scaling", v, &OPT_IMPLS)? { - let name = &v.ident; - etar.push(quote! { - Self::#name(eos) => eos.viscosity_reference(temperature, volume, moles) - }); - etac.push(quote! { - Self::#name(eos) => eos.viscosity_correlation(s_res, x) - }); - dr.push(quote! { - Self::#name(eos) => eos.diffusion_reference(temperature, volume, moles) - }); - dc.push(quote! { - Self::#name(eos) => eos.diffusion_correlation(s_res, x) - }); - thcr.push(quote! { - Self::#name(eos) => eos.thermal_conductivity_reference(temperature, volume, moles) - }); - thcc.push(quote! { - Self::#name(eos) => eos.thermal_conductivity_correlation(s_res, x) - }); - } - } - - Ok(quote! { - impl EntropyScaling for EosVariant { - fn viscosity_reference( - &self, - temperature: SINumber, - volume: SINumber, - moles: &SIArray1, - ) -> EosResult { - match self { - #(#etar,)* - _ => unimplemented!(), - } - } - - fn viscosity_correlation(&self, s_res: f64, x: &Array1) -> EosResult { - match self { - #(#etac,)* - _ => unimplemented!(), - } - } - - fn diffusion_reference( - &self, - temperature: SINumber, - volume: SINumber, - moles: &SIArray1, - ) -> EosResult { - match self { - #(#dr,)* - _ => unimplemented!(), - } - } - - fn diffusion_correlation(&self, s_res: f64, x: &Array1) -> EosResult { - match self { - #(#dc,)* - _ => unimplemented!(), - } - } - - fn thermal_conductivity_reference( - &self, - temperature: SINumber, - volume: SINumber, - moles: &SIArray1, - ) -> EosResult { - match self { - #(#thcr,)* - _ => unimplemented!(), - } - } - - fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1) -> EosResult { - match self { - #(#thcc,)* - _ => unimplemented!(), - } - } - } - }) -} diff --git a/feos-derive/src/functional_contribution.rs b/feos-derive/src/functional_contribution.rs index 0f6a25487..d000813e0 100644 --- a/feos-derive/src/functional_contribution.rs +++ b/feos-derive/src/functional_contribution.rs @@ -45,12 +45,12 @@ fn impl_functional_contribution( quote! { impl FunctionalContribution for #ident { - fn weight_functions + Copy+ScalarOperand>(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions + Copy+ScalarOperand>(&self, temperature: N) -> feos_dft::WeightFunctionInfo { match self { #(#weight_functions,)* } } - fn weight_functions_pdgt + Copy+ScalarOperand>(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy+ScalarOperand>(&self, temperature: N) -> feos_dft::WeightFunctionInfo { match self { #(#weight_functions_pdgt,)* } @@ -58,7 +58,7 @@ fn impl_functional_contribution( fn helmholtz_energy_density + Copy+ScalarOperand>( &self, temperature: N, - weighted_densities: ArrayView2, + weighted_densities: ndarray::ArrayView2, ) -> EosResult> { match self { #(#helmholtz_energy_density,)* diff --git a/feos-derive/src/lib.rs b/feos-derive/src/lib.rs index 1a0d586c1..2304446b4 100644 --- a/feos-derive/src/lib.rs +++ b/feos-derive/src/lib.rs @@ -16,6 +16,16 @@ mod functional_contribution; mod ideal_gas; mod residual; +// possible additional traits to implement +const OPT_IMPLS: [&str; 6] = [ + "molar_weight", + "entropy_scaling", + "functional", + "bond_lengths", + "fluid_parameters", + "pair_potential", +]; + fn implement(name: &str, variant: &syn::Variant, opts: &[&'static str]) -> syn::Result { let syn::Variant { attrs, .. } = variant; let mut implement = Ok(false); @@ -74,7 +84,7 @@ pub fn derive_residual(input: TokenStream) -> TokenStream { .into() } -#[proc_macro_derive(HelmholtzEnergyFunctional, attributes(implement))] +#[proc_macro_derive(HelmholtzEnergyFunctional, attributes(implement_dft))] pub fn derive_helmholtz_energy_functional(input: TokenStream) -> TokenStream { let input = parse_macro_input!(input as DeriveInput); expand_helmholtz_energy_functional(input) diff --git a/feos-derive/src/residual.rs b/feos-derive/src/residual.rs index 9d772be67..b70c5a47d 100644 --- a/feos-derive/src/residual.rs +++ b/feos-derive/src/residual.rs @@ -1,19 +1,16 @@ -use super::implement; +use super::{implement, OPT_IMPLS}; use quote::quote; use syn::DeriveInput; -// possible additional traits to implement -const OPT_IMPLS: [&str; 2] = ["molar_weight", "entropy_scaling"]; - pub(crate) fn expand_residual(input: DeriveInput) -> syn::Result { let variants = match input.data { syn::Data::Enum(syn::DataEnum { ref variants, .. }) => variants, _ => panic!("this derive macro only works on enums"), }; - let residual = impl_residual(variants); - let molar_weight = impl_molar_weight(variants)?; - let entropy_scaling = impl_entropy_scaling(variants)?; + let residual = impl_residual(&input.ident, variants); + let molar_weight = impl_molar_weight(&input.ident, variants)?; + let entropy_scaling = impl_entropy_scaling(&input.ident, variants)?; Ok(quote! { #residual #molar_weight @@ -22,6 +19,7 @@ pub(crate) fn expand_residual(input: DeriveInput) -> syn::Result, ) -> proc_macro2::TokenStream { let compute_max_density = variants.iter().map(|v| { @@ -38,7 +36,7 @@ fn impl_residual( }); quote! { - impl Residual for ResidualModel { + impl Residual for #ident { fn compute_max_density(&self, moles: &Array1) -> f64 { match self { #(#compute_max_density,)* @@ -54,38 +52,44 @@ fn impl_residual( } fn impl_molar_weight( + ident: &syn::Ident, variants: &syn::punctuated::Punctuated, ) -> syn::Result { let mut molar_weight = Vec::new(); let mut has_molar_weight = Vec::new(); for v in variants.iter() { + let name = &v.ident; if implement("molar_weight", v, &OPT_IMPLS)? { - let name = &v.ident; molar_weight.push(quote! { Self::#name(eos) => eos.molar_weight() }); has_molar_weight.push(quote! { Self::#name(_) => true }); + } else { + molar_weight.push(quote! { + Self::#name(eos) => panic!("{} does not provide molar weights and can not be used to calculate mass-specific properties", stringify!(#name)) + }); + has_molar_weight.push(quote! { + Self::#name(_) => false + }); } } Ok(quote! { - impl Molarweight for ResidualModel { + impl Molarweight for #ident { fn molar_weight(&self) -> MolarWeight> { match self { #(#molar_weight,)* - _ => unimplemented!(), } } } - impl ResidualModel { + impl #ident { pub fn has_molar_weight(&self) -> bool { match self { #(#has_molar_weight,)* - _ => false, } } } @@ -93,6 +97,7 @@ fn impl_molar_weight( } fn impl_entropy_scaling( + ident: &syn::Ident, variants: &syn::punctuated::Punctuated, ) -> syn::Result { let mut etar = Vec::new(); @@ -103,8 +108,8 @@ fn impl_entropy_scaling( let mut thcc = Vec::new(); for v in variants.iter() { + let name = &v.ident; if implement("entropy_scaling", v, &OPT_IMPLS)? { - let name = &v.ident; etar.push(quote! { Self::#name(eos) => eos.viscosity_reference(temperature, volume, moles) }); @@ -123,11 +128,30 @@ fn impl_entropy_scaling( thcc.push(quote! { Self::#name(eos) => eos.thermal_conductivity_correlation(s_res, x) }); + } else { + etar.push(quote! { + Self::#name(eos) => panic!("{} does not implement entropy scaling for transport properties!", stringify!(#name)) + }); + etac.push(quote! { + Self::#name(eos) => panic!("{} does not implement entropy scaling for transport properties!", stringify!(#name)) + }); + dr.push(quote! { + Self::#name(eos) => panic!("{} does not implement entropy scaling for transport properties!", stringify!(#name)) + }); + dc.push(quote! { + Self::#name(eos) => panic!("{} does not implement entropy scaling for transport properties!", stringify!(#name)) + }); + thcr.push(quote! { + Self::#name(eos) => panic!("{} does not implement entropy scaling for transport properties!", stringify!(#name)) + }); + thcc.push(quote! { + Self::#name(eos) => panic!("{} does not implement entropy scaling for transport properties!", stringify!(#name)) + }); } } Ok(quote! { - impl EntropyScaling for ResidualModel { + impl EntropyScaling for #ident { fn viscosity_reference( &self, temperature: Temperature, @@ -136,14 +160,12 @@ fn impl_entropy_scaling( ) -> EosResult { match self { #(#etar,)* - _ => unimplemented!(), } } fn viscosity_correlation(&self, s_res: f64, x: &Array1) -> EosResult { match self { #(#etac,)* - _ => unimplemented!(), } } @@ -155,14 +177,12 @@ fn impl_entropy_scaling( ) -> EosResult { match self { #(#dr,)* - _ => unimplemented!(), } } fn diffusion_correlation(&self, s_res: f64, x: &Array1) -> EosResult { match self { #(#dc,)* - _ => unimplemented!(), } } @@ -174,14 +194,12 @@ fn impl_entropy_scaling( ) -> EosResult { match self { #(#thcr,)* - _ => unimplemented!(), } } fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1) -> EosResult { match self { #(#thcc,)* - _ => unimplemented!(), } } } diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs index 447227fbc..a50f87f47 100644 --- a/feos-dft/src/adsorption/mod.rs +++ b/feos-dft/src/adsorption/mod.rs @@ -1,9 +1,9 @@ //! Adsorption profiles and isotherms. -use super::functional::{HelmholtzEnergyFunctional, DFT}; +use super::functional::HelmholtzEnergyFunctional; use super::solver::DFTSolver; use feos_core::{ - Components, Contributions, DensityInitialization, EosError, EosResult, ReferenceSystem, - Residual, SolverOptions, State, StateBuilder, + Contributions, DensityInitialization, EosError, EosResult, ReferenceSystem, SolverOptions, + State, StateBuilder, }; use ndarray::{Array1, Array2, Dimension, Ix1, Ix3, RemoveAxis}; use quantity::{Energy, MolarEnergy, Moles, Pressure, Temperature}; @@ -45,7 +45,7 @@ where D::Smaller: Dimension, ::Larger: Dimension, { - fn new(functional: &Arc>, profiles: Vec>>) -> Self { + fn new(functional: &Arc, profiles: Vec>>) -> Self { Self { components: functional.components(), profiles, @@ -54,7 +54,7 @@ where /// Calculate an adsorption isotherm (starting at low pressure) pub fn adsorption_isotherm>( - functional: &Arc>, + functional: &Arc, temperature: Temperature, pressure: &Pressure>, pore: &S, @@ -74,7 +74,7 @@ where /// Calculate an desorption isotherm (starting at high pressure) pub fn desorption_isotherm>( - functional: &Arc>, + functional: &Arc, temperature: Temperature, pressure: &Pressure>, pore: &S, @@ -99,7 +99,7 @@ where /// Calculate an equilibrium isotherm pub fn equilibrium_isotherm>( - functional: &Arc>, + functional: &Arc, temperature: Temperature, pressure: &Pressure>, pore: &S, @@ -182,7 +182,7 @@ where } fn isotherm>( - functional: &Arc>, + functional: &Arc, temperature: Temperature, pressure: &Pressure>, pore: &S, @@ -245,7 +245,7 @@ where /// Calculate the phase transition from an empty to a filled pore. #[expect(clippy::too_many_arguments)] pub fn phase_equilibrium>( - functional: &Arc>, + functional: &Arc, temperature: Temperature, p_min: Pressure, p_max: Pressure, diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index b90e91fad..c0b768308 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -1,12 +1,14 @@ use crate::adsorption::{ExternalPotential, FluidParameters}; use crate::convolver::ConvolverFFT; -use crate::functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; +use crate::functional::{HelmholtzEnergyFunctional, MoleculeShape}; use crate::functional_contribution::FunctionalContribution; use crate::geometry::{Axis, Geometry, Grid}; use crate::profile::{DFTProfile, MAX_POTENTIAL}; use crate::solver::DFTSolver; use crate::WeightFunctionInfo; -use feos_core::{Components, Contributions, EosResult, ReferenceSystem, State, StateBuilder}; +use feos_core::{ + Components, Contributions, EosResult, ReferenceSystem, Residual, State, StateBuilder, StateHD, +}; use ndarray::{prelude::*, ScalarOperand}; use ndarray::{Axis as Axis_nd, RemoveAxis}; use num_dual::linalg::LU; @@ -58,7 +60,7 @@ pub trait PoreSpecification { /// Initialize a new single pore. fn initialize( &self, - bulk: &State>, + bulk: &State, density: Option<&Density>>, external_potential: Option<&Array>, ) -> EosResult>; @@ -129,7 +131,7 @@ where Ok(self) } - pub fn update_bulk(mut self, bulk: &State>) -> Self { + pub fn update_bulk(mut self, bulk: &State) -> Self { self.profile.bulk = bulk.clone(); self.grand_potential = None; self.interfacial_tension = None; @@ -194,7 +196,7 @@ where impl PoreSpecification for Pore1D { fn initialize( &self, - bulk: &State>, + bulk: &State, density: Option<&Density>>, external_potential: Option<&Array2>, ) -> EosResult> { @@ -308,10 +310,10 @@ struct Helium { } impl Helium { - fn new() -> DFT { + fn new() -> Self { let epsilon = arr1(&[EPSILON_HE]); let sigma = arr1(&[SIGMA_HE]); - DFT(Self { epsilon, sigma }) + Self { epsilon, sigma } } } @@ -325,6 +327,19 @@ impl Components for Helium { } } +impl Residual for Helium { + fn compute_max_density(&self, _: &Array1) -> f64 { + 1.0 + } + + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + self.evaluate_bulk(state) + } +} + impl HelmholtzEnergyFunctional for Helium { type Contribution = HeliumContribution; @@ -332,10 +347,6 @@ impl HelmholtzEnergyFunctional for Helium { Box::new([].into_iter()) } - fn compute_max_density(&self, _: &Array1) -> f64 { - 1.0 - } - fn molecule_shape(&self) -> MoleculeShape { MoleculeShape::Spherical(1) } diff --git a/feos-dft/src/adsorption/pore2d.rs b/feos-dft/src/adsorption/pore2d.rs index af5614e44..6b95d7e29 100644 --- a/feos-dft/src/adsorption/pore2d.rs +++ b/feos-dft/src/adsorption/pore2d.rs @@ -1,5 +1,5 @@ use super::{FluidParameters, PoreProfile, PoreSpecification}; -use crate::{Axis, DFTProfile, Grid, HelmholtzEnergyFunctional, DFT}; +use crate::{Axis, DFTProfile, Grid, HelmholtzEnergyFunctional}; use feos_core::{EosResult, State}; use ndarray::{Array3, Ix2}; use quantity::{Angle, Density, Length}; @@ -25,7 +25,7 @@ impl Pore2D { impl PoreSpecification for Pore2D { fn initialize( &self, - bulk: &State>, + bulk: &State, density: Option<&Density>>, external_potential: Option<&Array3>, ) -> EosResult> { diff --git a/feos-dft/src/adsorption/pore3d.rs b/feos-dft/src/adsorption/pore3d.rs index 02a1238e8..e03800521 100644 --- a/feos-dft/src/adsorption/pore3d.rs +++ b/feos-dft/src/adsorption/pore3d.rs @@ -1,6 +1,6 @@ use super::pore::{PoreProfile, PoreSpecification}; use crate::adsorption::FluidParameters; -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::{Axis, Grid}; use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL}; use feos_core::{EosError, EosResult, ReferenceSystem, State}; @@ -51,7 +51,7 @@ pub type PoreProfile3D = PoreProfile; impl PoreSpecification for Pore3D { fn initialize( &self, - bulk: &State>, + bulk: &State, density: Option<&Density>>, external_potential: Option<&Array4>, ) -> EosResult> { diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 4f8460303..db041a957 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -4,18 +4,17 @@ use crate::functional_contribution::*; use crate::ideal_chain_contribution::IdealChainContribution; use crate::solvation::PairPotential; use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; -use feos_core::{Components, EosResult, EquationOfState, IdealGas, Molarweight, Residual, StateHD}; +use feos_core::{EosResult, EquationOfState, IdealGas, Residual, StateHD}; use ndarray::*; use num_dual::*; use petgraph::graph::{Graph, UnGraph}; use petgraph::visit::EdgeRef; use petgraph::Directed; -use quantity::MolarWeight; use std::borrow::Cow; -use std::ops::{Deref, MulAssign}; +use std::ops::MulAssign; use std::sync::Arc; -impl HelmholtzEnergyFunctional +impl HelmholtzEnergyFunctional for EquationOfState { type Contribution = F::Contribution; @@ -28,10 +27,6 @@ impl HelmholtzEnergyF self.residual.molecule_shape() } - fn compute_max_density(&self, moles: &Array1) -> f64 { - self.residual.compute_max_density(moles) - } - fn bond_lengths + Copy>(&self, temperature: N) -> UnGraph<(), N> { self.residual.bond_lengths(temperature) } @@ -43,7 +38,7 @@ impl PairPotential for EquationOfState { } } -impl FluidParameters for EquationOfState { +impl FluidParameters for EquationOfState { fn epsilon_k_ff(&self) -> Array1 { self.residual.epsilon_k_ff() } @@ -53,80 +48,6 @@ impl FluidParameters for Equati } } -/// Wrapper struct for the [HelmholtzEnergyFunctional] trait. -/// -/// Needed (for now) to generically implement the `Residual` -/// trait for Helmholtz energy functionals. -#[derive(Clone)] -pub struct DFT(pub F); - -impl DFT { - pub fn into>(self) -> DFT { - DFT(self.0.into()) - } -} - -impl Deref for DFT { - type Target = F; - fn deref(&self) -> &F { - &self.0 - } -} - -impl DFT { - pub fn ideal_gas(self, ideal_gas: I) -> DFT> { - DFT(EquationOfState::new(Arc::new(ideal_gas), Arc::new(self.0))) - } -} - -impl Components for DFT { - fn components(&self) -> usize { - self.0.components() - } - - fn subset(&self, component_list: &[usize]) -> Self { - Self(self.0.subset(component_list)) - } -} - -impl Residual for DFT { - fn compute_max_density(&self, moles: &Array1) -> f64 { - self.0.compute_max_density(moles) - } - - fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( - &self, - state: &StateHD, - ) -> Vec<(String, D)> { - let mut res: Vec<(String, D)> = self - .0 - .contributions() - .map(|c| (c.to_string(), c.helmholtz_energy(state))) - .collect(); - res.push(( - self.ideal_chain_contribution().to_string(), - self.ideal_chain_contribution().helmholtz_energy(state), - )); - res - } -} - -impl Molarweight for DFT { - fn molar_weight(&self) -> MolarWeight> { - self.0.molar_weight() - } -} - -impl IdealGas for DFT { - fn ln_lambda3 + Copy>(&self, temperature: D) -> Array1 { - self.0.ln_lambda3(temperature) - } - - fn ideal_gas_model(&self) -> String { - self.0.ideal_gas_model() - } -} - /// Different representations for molecules within DFT. pub enum MoleculeShape<'a> { /// For spherical molecules, the number of components. @@ -140,7 +61,7 @@ pub enum MoleculeShape<'a> { } /// A general Helmholtz energy functional. -pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { +pub trait HelmholtzEnergyFunctional: Residual + Sized { type Contribution: FunctionalContribution; /// Return a slice of [FunctionalContribution]s. @@ -149,14 +70,6 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { /// Return the shape of the molecules and the necessary specifications. fn molecule_shape(&self) -> MoleculeShape; - /// Return the maximum density in Angstrom^-3. - /// - /// This value is used as an estimate for a liquid phase for phase - /// equilibria and other iterations. It is not explicitly meant to - /// be a mathematical limit for the density (if those exist in the - /// equation of state anyways). - fn compute_max_density(&self, moles: &Array1) -> f64; - /// Overwrite this, if the functional consists of heterosegmented chains. fn bond_lengths + Copy>(&self, _temperature: N) -> UnGraph<(), N> { Graph::with_capacity(0, 0) @@ -307,4 +220,19 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { i } + + fn evaluate_bulk + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + let mut res: Vec<(String, D)> = self + .contributions() + .map(|c| (c.to_string(), c.helmholtz_energy(state))) + .collect(); + res.push(( + self.ideal_chain_contribution().to_string(), + self.ideal_chain_contribution().helmholtz_energy(state), + )); + res + } } diff --git a/feos-dft/src/interface/mod.rs b/feos-dft/src/interface/mod.rs index 39d7c81ee..c9ab3cee2 100644 --- a/feos-dft/src/interface/mod.rs +++ b/feos-dft/src/interface/mod.rs @@ -1,6 +1,7 @@ //! Density profiles at planar interfaces and interfacial tensions. -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::{Axis, Grid}; +use crate::pdgt::PdgtFunctionalProperties; use crate::profile::{DFTProfile, DFTSpecifications}; use crate::solver::DFTSolver; use feos_core::{Contributions, EosError, EosResult, PhaseEquilibrium, ReferenceSystem}; @@ -16,7 +17,7 @@ const MIN_WIDTH: f64 = 100.0; /// Density profile and properties of a planar interface. pub struct PlanarInterface { pub profile: DFTProfile, - pub vle: PhaseEquilibrium, 2>, + pub vle: PhaseEquilibrium, pub surface_tension: Option, pub equimolar_radius: Option, } @@ -62,7 +63,7 @@ impl PlanarInterface { } impl PlanarInterface { - pub fn new(vle: &PhaseEquilibrium, 2>, n_grid: usize, l_grid: Length) -> Self { + pub fn new(vle: &PhaseEquilibrium, n_grid: usize, l_grid: Length) -> Self { // generate grid let grid = Grid::Cartesian1(Axis::new_cartesian(n_grid, l_grid, None)); @@ -75,7 +76,7 @@ impl PlanarInterface { } pub fn from_tanh( - vle: &PhaseEquilibrium, 2>, + vle: &PhaseEquilibrium, n_grid: usize, l_grid: Length, critical_temperature: Temperature, @@ -111,7 +112,7 @@ impl PlanarInterface { } pub fn from_pdgt( - vle: &PhaseEquilibrium, 2>, + vle: &PhaseEquilibrium, n_grid: usize, fix_equimolar_surface: bool, ) -> EosResult { @@ -341,10 +342,10 @@ impl PlanarInterface { } fn interp_symmetric( - vle_pdgt: &PhaseEquilibrium, 2>, + vle_pdgt: &PhaseEquilibrium, z_pdgt: Length>, rho_pdgt: Density>, - vle: &PhaseEquilibrium, 2>, + vle: &PhaseEquilibrium, z: &Array1, radius: Length, ) -> EosResult>> { diff --git a/feos-dft/src/interface/surface_tension_diagram.rs b/feos-dft/src/interface/surface_tension_diagram.rs index 9f4db6ef3..65741c751 100644 --- a/feos-dft/src/interface/surface_tension_diagram.rs +++ b/feos-dft/src/interface/surface_tension_diagram.rs @@ -1,5 +1,5 @@ use super::PlanarInterface; -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::HelmholtzEnergyFunctional; use crate::solver::DFTSolver; use feos_core::{PhaseEquilibrium, ReferenceSystem, StateVec}; use ndarray::{Array1, Array2}; @@ -15,7 +15,7 @@ pub struct SurfaceTensionDiagram { // #[expect(clippy::ptr_arg)] impl SurfaceTensionDiagram { pub fn new( - dia: &[PhaseEquilibrium, 2>], + dia: &[PhaseEquilibrium], init_densities: Option, n_grid: Option, l_grid: Option, @@ -67,11 +67,11 @@ impl SurfaceTensionDiagram { Self { profiles } } - pub fn vapor(&self) -> StateVec<'_, DFT> { + pub fn vapor(&self) -> StateVec<'_, F> { self.profiles.iter().map(|p| p.vle.vapor()).collect() } - pub fn liquid(&self) -> StateVec<'_, DFT> { + pub fn liquid(&self) -> StateVec<'_, F> { self.profiles.iter().map(|p| p.vle.liquid()).collect() } diff --git a/feos-dft/src/lib.rs b/feos-dft/src/lib.rs index c2d78fcc3..ee08dee8d 100644 --- a/feos-dft/src/lib.rs +++ b/feos-dft/src/lib.rs @@ -15,9 +15,10 @@ mod solver; mod weight_functions; pub use convolver::{Convolver, ConvolverFFT}; -pub use functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; +pub use functional::{HelmholtzEnergyFunctional, MoleculeShape}; pub use functional_contribution::FunctionalContribution; pub use geometry::{Axis, Geometry, Grid}; +pub use pdgt::PdgtFunctionalProperties; pub use profile::{DFTProfile, DFTSpecification, DFTSpecifications}; pub use solver::{DFTSolver, DFTSolverLog}; pub use weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; diff --git a/feos-dft/src/pdgt.rs b/feos-dft/src/pdgt.rs index 83e423e64..7692b06de 100644 --- a/feos-dft/src/pdgt.rs +++ b/feos-dft/src/pdgt.rs @@ -1,7 +1,7 @@ -use super::functional::{HelmholtzEnergyFunctional, DFT}; +use super::functional::HelmholtzEnergyFunctional; use super::functional_contribution::FunctionalContribution; use super::weight_functions::WeightFunctionInfo; -use feos_core::{Components, Contributions, EosResult, PhaseEquilibrium, ReferenceSystem}; +use feos_core::{Contributions, EosResult, PhaseEquilibrium, ReferenceSystem}; use ndarray::*; use num_dual::Dual2_64; use quantity::{ @@ -144,8 +144,9 @@ trait PdgtProperties: FunctionalContribution { impl PdgtProperties for T {} -impl DFT { - pub fn solve_pdgt( +pub trait PdgtFunctionalProperties: HelmholtzEnergyFunctional { + // impl T { + fn solve_pdgt( &self, vle: &PhaseEquilibrium, n_grid: usize, @@ -244,6 +245,8 @@ impl DFT { } } +impl PdgtFunctionalProperties for T {} + fn gradient, UX: Copy>( df: &Quantity, UF>, dx: Quantity, diff --git a/feos-dft/src/profile/mod.rs b/feos-dft/src/profile/mod.rs index 230ff45ce..0c7acf5b0 100644 --- a/feos-dft/src/profile/mod.rs +++ b/feos-dft/src/profile/mod.rs @@ -1,8 +1,8 @@ use crate::convolver::{BulkConvolver, Convolver, ConvolverFFT}; -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::Grid; use crate::solver::{DFTSolver, DFTSolverLog}; -use feos_core::{Components, EosError, EosResult, ReferenceSystem, State}; +use feos_core::{EosError, EosResult, ReferenceSystem, State}; use ndarray::{ Array, Array1, Array2, Array3, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3, RemoveAxis, @@ -101,12 +101,12 @@ impl DFTSpecification for DFTS pub struct DFTProfile { pub grid: Grid, pub convolver: Arc>, - pub dft: Arc>, + pub dft: Arc, pub temperature: Temperature, pub density: Density>, pub specification: Arc>, pub external_potential: Array, - pub bulk: State>, + pub bulk: State, pub solver_log: Option, pub lanczos: Option, } @@ -205,7 +205,7 @@ where /// after this call if something else is required. pub fn new( grid: Grid, - bulk: &State>, + bulk: &State, external_potential: Option>, density: Option<&Density>>, lanczos: Option, diff --git a/feos-dft/src/solvation/pair_correlation.rs b/feos-dft/src/solvation/pair_correlation.rs index 2d5d2400d..4320976ed 100644 --- a/feos-dft/src/solvation/pair_correlation.rs +++ b/feos-dft/src/solvation/pair_correlation.rs @@ -1,5 +1,5 @@ //! Functionalities for the calculation of pair correlation functions. -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::HelmholtzEnergyFunctional; use crate::profile::MAX_POTENTIAL; use crate::solver::DFTSolver; use crate::{Axis, DFTProfile, Grid}; @@ -34,7 +34,7 @@ impl Clone for PairCorrelation { } impl PairCorrelation { - pub fn new(bulk: &State>, test_particle: usize, n_grid: usize, width: Length) -> Self { + pub fn new(bulk: &State, test_particle: usize, n_grid: usize, width: Length) -> Self { let dft = &bulk.eos; // generate grid diff --git a/feos-dft/src/solvation/solvation_profile.rs b/feos-dft/src/solvation/solvation_profile.rs index eff71bc27..af3e8263c 100644 --- a/feos-dft/src/solvation/solvation_profile.rs +++ b/feos-dft/src/solvation/solvation_profile.rs @@ -1,5 +1,5 @@ use crate::adsorption::FluidParameters; -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::{Axis, Grid}; use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL}; use crate::solver::DFTSolver; @@ -52,7 +52,7 @@ impl SolvationProfile { impl SolvationProfile { #[expect(clippy::too_many_arguments)] pub fn new( - bulk: &State>, + bulk: &State, n_grid: [usize; 3], coordinates: Length>, sigma_ss: Array1, diff --git a/src/eos.rs b/src/eos.rs index cd795cb46..e13ad5a17 100644 --- a/src/eos.rs +++ b/src/eos.rs @@ -1,22 +1,10 @@ -#[cfg(feature = "epcsaft")] -use crate::epcsaft::ElectrolytePcSaft; -#[cfg(feature = "gc_pcsaft")] -use crate::gc_pcsaft::GcPcSaft; -#[cfg(feature = "pcsaft")] -use crate::pcsaft::PcSaft; -#[cfg(feature = "pets")] -use crate::pets::Pets; -#[cfg(feature = "saftvrmie")] -use crate::saftvrmie::SaftVRMie; -#[cfg(feature = "saftvrqmie")] -use crate::saftvrqmie::SaftVRQMie; -#[cfg(feature = "uvtheory")] -use crate::uvtheory::UVTheory; use feos_core::cubic::PengRobinson; -#[cfg(feature = "python")] -use feos_core::python::user_defined::PyResidual; use feos_core::*; use feos_derive::{Components, Residual}; +#[cfg(feature = "dft")] +use feos_derive::{FunctionalContribution, HelmholtzEnergyFunctional}; +#[cfg(feature = "dft")] +use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional}; use ndarray::{Array1, ScalarOperand}; use num_dual::DualNum; use quantity::*; @@ -25,33 +13,77 @@ use quantity::*; /// /// Particularly relevant for situations in which generic types /// are undesirable (e.g. FFI). +#[cfg_attr(feature = "dft", derive(HelmholtzEnergyFunctional))] #[derive(Components, Residual)] pub enum ResidualModel { + // Equations of state NoResidual(NoResidual), #[cfg(feature = "pcsaft")] #[implement(entropy_scaling, molar_weight)] - PcSaft(PcSaft), + PcSaft(crate::pcsaft::PcSaft), + #[cfg(feature = "epcsaft")] #[implement(molar_weight)] - ElectrolytePcSaft(ElectrolytePcSaft), + ElectrolytePcSaft(crate::epcsaft::ElectrolytePcSaft), + #[cfg(feature = "gc_pcsaft")] #[implement(molar_weight)] - GcPcSaft(GcPcSaft), + GcPcSaft(crate::gc_pcsaft::GcPcSaft), + #[implement(molar_weight)] PengRobinson(PengRobinson), + #[cfg(feature = "python")] #[implement(molar_weight)] - Python(PyResidual), + Python(feos_core::python::user_defined::PyResidual), + #[cfg(feature = "saftvrqmie")] #[implement(entropy_scaling, molar_weight)] - SaftVRQMie(SaftVRQMie), + SaftVRQMie(crate::saftvrqmie::SaftVRQMie), + #[cfg(feature = "saftvrmie")] #[implement(molar_weight)] - SaftVRMie(SaftVRMie), + SaftVRMie(crate::saftvrmie::SaftVRMie), + #[cfg(feature = "pets")] #[implement(molar_weight)] - Pets(Pets), + Pets(crate::pets::Pets), + #[cfg(feature = "uvtheory")] #[implement(molar_weight)] - UVTheory(UVTheory), + UVTheory(crate::uvtheory::UVTheory), + + // Helmholtz energy functionals + #[cfg(feature = "dft")] + #[implement(molar_weight, functional, fluid_parameters, pair_potential)] + PcSaftFunctional(crate::pcsaft::PcSaftFunctional), + + #[cfg(feature = "gc_pcsaft")] + #[implement(molar_weight, functional, fluid_parameters, bond_lengths)] + GcPcSaftFunctional(crate::gc_pcsaft::GcPcSaftFunctional), + + #[cfg(feature = "pets")] + #[implement(molar_weight, functional, fluid_parameters, pair_potential)] + PetsFunctional(crate::pets::PetsFunctional), + + #[implement(functional, fluid_parameters, pair_potential)] + FmtFunctional(crate::hard_sphere::FMTFunctional), + + #[cfg(feature = "saftvrqmie")] + #[implement(molar_weight, functional, fluid_parameters, pair_potential)] + SaftVRQMieFunctional(crate::saftvrqmie::SaftVRQMieFunctional), +} + +#[cfg(feature = "dft")] +#[derive(FunctionalContribution)] +pub enum FunctionalContributionVariant { + #[cfg(feature = "pcsaft")] + PcSaftFunctional(crate::pcsaft::PcSaftFunctionalContribution), + #[cfg(feature = "gc_pcsaft")] + GcPcSaftFunctional(crate::gc_pcsaft::GcPcSaftFunctionalContribution), + #[cfg(feature = "pets")] + PetsFunctional(crate::pets::PetsFunctionalContribution), + Fmt(crate::hard_sphere::FMTContribution), + #[cfg(feature = "saftvrqmie")] + SaftVRQMieFunctional(crate::saftvrqmie::SaftVRQMieFunctionalContribution), } diff --git a/src/functional.rs b/src/functional.rs deleted file mode 100644 index 784ba4166..000000000 --- a/src/functional.rs +++ /dev/null @@ -1,55 +0,0 @@ -#[cfg(feature = "gc_pcsaft")] -use crate::gc_pcsaft::{GcPcSaftFunctional, GcPcSaftFunctionalContribution}; -use crate::hard_sphere::{FMTContribution, FMTFunctional, HardSphereParameters}; -#[cfg(feature = "pcsaft")] -use crate::pcsaft::{PcSaftFunctional, PcSaftFunctionalContribution}; -#[cfg(feature = "pets")] -use crate::pets::{PetsFunctional, PetsFunctionalContribution}; -#[cfg(feature = "saftvrqmie")] -use crate::saftvrqmie::{SaftVRQMieFunctional, SaftVRQMieFunctionalContribution}; -use quantity::MolarWeight; -use feos_core::*; -use feos_derive::FunctionalContribution; -use feos_derive::{Components, HelmholtzEnergyFunctional}; -use feos_dft::adsorption::*; -use feos_dft::solvation::*; -use feos_dft::*; -use ndarray::{Array1, Array2, ArrayView2, ScalarOperand}; -use num_dual::DualNum; -use petgraph::graph::UnGraph; -use petgraph::Graph; - -/// Collection of different [HelmholtzEnergyFunctional] implementations. -/// -/// Particularly relevant for situations in which generic types -/// are undesirable (e.g. FFI). -#[derive(Components, HelmholtzEnergyFunctional)] -pub enum FunctionalVariant { - #[cfg(feature = "pcsaft")] - #[implement(fluid_parameters, molar_weight, pair_potential)] - PcSaft(PcSaftFunctional), - #[cfg(feature = "gc_pcsaft")] - #[implement(fluid_parameters, molar_weight, bond_lengths)] - GcPcSaft(GcPcSaftFunctional), - #[cfg(feature = "pets")] - #[implement(fluid_parameters, molar_weight, pair_potential)] - Pets(PetsFunctional), - #[implement(fluid_parameters, pair_potential)] - Fmt(FMTFunctional), - #[cfg(feature = "saftvrqmie")] - #[implement(fluid_parameters, molar_weight, pair_potential)] - SaftVRQMie(SaftVRQMieFunctional), -} - -#[derive(FunctionalContribution)] -pub enum FunctionalContributionVariant { - #[cfg(feature = "pcsaft")] - PcSaftFunctional(PcSaftFunctionalContribution), - #[cfg(feature = "gc_pcsaft")] - GcPcSaftFunctional(GcPcSaftFunctionalContribution), - #[cfg(feature = "pets")] - PetsFunctional(PetsFunctionalContribution), - Fmt(FMTContribution), - #[cfg(feature = "saftvrqmie")] - SaftVRQMieFunctional(SaftVRQMieFunctionalContribution), -} diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 50191c9cd..a4bf65fbf 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -3,13 +3,11 @@ use super::record::GcPcSaftAssociationRecord; use crate::association::{Association, AssociationStrength}; use crate::hard_sphere::{FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; use feos_core::parameter::ParameterHetero; -use feos_core::{Components, EosResult, Molarweight}; +use feos_core::{Components, EosResult, Molarweight, Residual, StateHD}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; -use feos_dft::{ - FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, -}; -use ndarray::{Array1, ArrayView2, ScalarOperand}; +use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape}; +use ndarray::{Array1, ScalarOperand}; use num_dual::DualNum; use petgraph::graph::UnGraph; use quantity::{MolarWeight, GRAM, MOL}; @@ -31,7 +29,7 @@ pub struct GcPcSaftFunctional { } impl GcPcSaftFunctional { - pub fn new(parameters: Arc) -> DFT { + pub fn new(parameters: Arc) -> Self { Self::with_options( parameters, FMTVersion::WhiteBear, @@ -43,12 +41,12 @@ impl GcPcSaftFunctional { parameters: Arc, fmt_version: FMTVersion, saft_options: GcPcSaftOptions, - ) -> DFT { - DFT(Self { + ) -> Self { + Self { parameters, fmt_version, options: saft_options, - }) + } } } @@ -63,17 +61,10 @@ impl Components for GcPcSaftFunctional { self.fmt_version, self.options, ) - .0 } } -impl HelmholtzEnergyFunctional for GcPcSaftFunctional { - type Contribution = GcPcSaftFunctionalContribution; - - fn molecule_shape(&self) -> MoleculeShape { - MoleculeShape::Heterosegmented(&self.parameters.component_index) - } - +impl Residual for GcPcSaftFunctional { fn compute_max_density(&self, moles: &Array1) -> f64 { let p = &self.parameters; let moles_segments: Array1 = p.component_index.iter().map(|&i| moles[i]).collect(); @@ -81,6 +72,21 @@ impl HelmholtzEnergyFunctional for GcPcSaftFunctional { / (FRAC_PI_6 * &p.m * p.sigma.mapv(|v| v.powi(3)) * moles_segments).sum() } + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + self.evaluate_bulk(state) + } +} + +impl HelmholtzEnergyFunctional for GcPcSaftFunctional { + type Contribution = GcPcSaftFunctionalContribution; + + fn molecule_shape(&self) -> MoleculeShape { + MoleculeShape::Heterosegmented(&self.parameters.component_index) + } + fn contributions(&self) -> Box> { let mut contributions = Vec::with_capacity(4); diff --git a/src/gc_pcsaft/micelles.rs b/src/gc_pcsaft/micelles.rs index 671b7c35e..0033865f7 100644 --- a/src/gc_pcsaft/micelles.rs +++ b/src/gc_pcsaft/micelles.rs @@ -132,7 +132,7 @@ impl MicelleProfile { impl MicelleProfile { fn new( - bulk: &State>, + bulk: &State, axis: Axis, initialization: MicelleInitialization, specification: MicelleSpecification, @@ -185,7 +185,7 @@ impl MicelleProfile { } pub fn new_spherical( - bulk: &State>, + bulk: &State, n_grid: usize, width: SINumber, initialization: MicelleInitialization, @@ -200,7 +200,7 @@ impl MicelleProfile { } pub fn new_cylindrical( - bulk: &State>, + bulk: &State, n_grid: usize, width: SINumber, initialization: MicelleInitialization, diff --git a/src/hard_sphere/dft.rs b/src/hard_sphere/dft.rs index 8be7677a9..0d053f6a2 100644 --- a/src/hard_sphere/dft.rs +++ b/src/hard_sphere/dft.rs @@ -1,9 +1,9 @@ -use feos_core::{Components, EosResult}; +use feos_core::{Components, EosResult, Residual, StateHD}; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; use feos_dft::{ FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunction, - WeightFunctionInfo, WeightFunctionShape, DFT, + WeightFunctionInfo, WeightFunctionShape, }; use ndarray::*; use num_dual::DualNum; @@ -313,14 +313,14 @@ pub struct FMTFunctional { } impl FMTFunctional { - pub fn new(sigma: &Array1, version: FMTVersion) -> DFT { + pub fn new(sigma: &Array1, version: FMTVersion) -> Self { let properties = Arc::new(HardSphereParameters { sigma: sigma.clone(), }); - DFT(Self { + Self { properties, version, - }) + } } } @@ -334,7 +334,20 @@ impl Components for FMTFunctional { .iter() .map(|&c| self.properties.sigma[c]) .collect(); - Self::new(&sigma, self.version).0 + Self::new(&sigma, self.version) + } +} + +impl Residual for FMTFunctional { + fn compute_max_density(&self, moles: &Array1) -> f64 { + moles.sum() / (moles * &self.properties.sigma).sum() * 1.2 + } + + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + self.evaluate_bulk(state) } } @@ -348,10 +361,6 @@ impl HelmholtzEnergyFunctional for FMTFunctional { ))) } - 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()) } diff --git a/src/lib.rs b/src/lib.rs index 39f741a3f..fdc944708 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -35,10 +35,6 @@ #![warn(clippy::all)] #![warn(clippy::allow_attributes)] -#[cfg(feature = "dft")] -mod functional; -#[cfg(feature = "dft")] -pub use functional::FunctionalVariant; mod eos; pub use eos::ResidualModel; diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index d2d186960..f7d95657d 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -3,14 +3,12 @@ use crate::association::Association; use crate::hard_sphere::{FMTContribution, FMTVersion}; use crate::pcsaft::eos::PcSaftOptions; use feos_core::parameter::Parameter; -use feos_core::{Components, EosResult, Molarweight}; +use feos_core::{Components, EosResult, Molarweight, Residual, StateHD}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; -use feos_dft::{ - FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, -}; -use ndarray::{Array1, Array2, ArrayView2, ScalarOperand}; +use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape}; +use ndarray::{Array1, Array2, ScalarOperand}; use num_dual::DualNum; use num_traits::One; use quantity::{MolarWeight, GRAM, MOL}; @@ -33,11 +31,11 @@ pub struct PcSaftFunctional { } impl PcSaftFunctional { - pub fn new(parameters: Arc) -> DFT { + pub fn new(parameters: Arc) -> Self { Self::with_options(parameters, FMTVersion::WhiteBear, PcSaftOptions::default()) } - pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> DFT { + pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> Self { Self::with_options(parameters, fmt_version, PcSaftOptions::default()) } @@ -45,12 +43,12 @@ impl PcSaftFunctional { parameters: Arc, fmt_version: FMTVersion, saft_options: PcSaftOptions, - ) -> DFT { - DFT(Self { + ) -> Self { + Self { parameters, fmt_version, options: saft_options, - }) + } } } @@ -65,19 +63,27 @@ impl Components for PcSaftFunctional { self.fmt_version, self.options, ) - .0 } } -impl HelmholtzEnergyFunctional for PcSaftFunctional { - type Contribution = PcSaftFunctionalContribution; - +impl Residual for PcSaftFunctional { fn compute_max_density(&self, moles: &Array1) -> f64 { self.options.max_eta * moles.sum() / (FRAC_PI_6 * &self.parameters.m * self.parameters.sigma.mapv(|v| v.powi(3)) * moles) .sum() } + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + self.evaluate_bulk(state) + } +} + +impl HelmholtzEnergyFunctional for PcSaftFunctional { + type Contribution = PcSaftFunctionalContribution; + fn contributions(&self) -> Box> { let mut contributions = Vec::with_capacity(4); diff --git a/src/pets/dft/mod.rs b/src/pets/dft/mod.rs index bc76ae905..8c8b0920c 100644 --- a/src/pets/dft/mod.rs +++ b/src/pets/dft/mod.rs @@ -3,14 +3,12 @@ use super::parameters::PetsParameters; use crate::hard_sphere::{FMTContribution, FMTVersion}; use dispersion::AttractiveFunctional; use feos_core::parameter::Parameter; -use feos_core::{Components, EosResult, Molarweight}; +use feos_core::{Components, EosResult, Molarweight, Residual, StateHD}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; -use feos_dft::{ - FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, -}; -use ndarray::{Array1, Array2, ArrayView2, ScalarOperand}; +use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape}; +use ndarray::{Array1, Array2, ScalarOperand}; use num_dual::DualNum; use pure_pets_functional::*; use quantity::{MolarWeight, GRAM, MOL}; @@ -33,12 +31,12 @@ impl PetsFunctional { /// /// # Defaults /// `FMTVersion`: `FMTVersion::WhiteBear` - pub fn new(parameters: Arc) -> DFT { + pub fn new(parameters: Arc) -> Self { Self::with_options(parameters, FMTVersion::WhiteBear, PetsOptions::default()) } /// PeTS functional with default options for and provided FMT version. - pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> DFT { + pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> Self { Self::with_options(parameters, fmt_version, PetsOptions::default()) } @@ -47,12 +45,12 @@ impl PetsFunctional { parameters: Arc, fmt_version: FMTVersion, pets_options: PetsOptions, - ) -> DFT { - DFT(Self { + ) -> Self { + Self { parameters, fmt_version, options: pets_options, - }) + } } } @@ -67,7 +65,20 @@ impl Components for PetsFunctional { self.fmt_version, self.options, ) - .0 + } +} + +impl Residual for PetsFunctional { + fn compute_max_density(&self, moles: &Array1) -> f64 { + self.options.max_eta * moles.sum() + / (FRAC_PI_6 * self.parameters.sigma.mapv(|v| v.powi(3)) * moles).sum() + } + + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + self.evaluate_bulk(state) } } @@ -78,11 +89,6 @@ impl HelmholtzEnergyFunctional for PetsFunctional { MoleculeShape::Spherical(self.parameters.sigma.len()) } - fn compute_max_density(&self, moles: &Array1) -> f64 { - self.options.max_eta * moles.sum() - / (FRAC_PI_6 * self.parameters.sigma.mapv(|v| v.powi(3)) * moles).sum() - } - fn contributions(&self) -> Box<(dyn Iterator)> { let mut contributions = Vec::with_capacity(2); diff --git a/src/python/dft.rs b/src/python/dft.rs index 6536a8800..2d03735bf 100644 --- a/src/python/dft.rs +++ b/src/python/dft.rs @@ -1,14 +1,9 @@ -#[cfg(feature = "estimator")] -use crate::estimator::*; -use crate::functional::FunctionalVariant; #[cfg(feature = "gc_pcsaft")] use crate::gc_pcsaft::python::PyGcPcSaftFunctionalParameters; #[cfg(feature = "gc_pcsaft")] use crate::gc_pcsaft::{GcPcSaftFunctional, GcPcSaftOptions}; use crate::hard_sphere::{FMTFunctional, FMTVersion}; use crate::ideal_gas::IdealGasModel; -#[cfg(feature = "estimator")] -use crate::impl_estimator; #[cfg(feature = "pcsaft")] use crate::pcsaft::python::PyPcSaftParameters; #[cfg(feature = "pcsaft")] @@ -21,7 +16,9 @@ use crate::pets::{PetsFunctional, PetsOptions}; use crate::saftvrqmie::python::PySaftVRQMieParameters; #[cfg(feature = "saftvrqmie")] use crate::saftvrqmie::{SaftVRQMieFunctional, SaftVRQMieOptions}; +use crate::ResidualModel; +use super::eos::{PyEquationOfState, PyPhaseEquilibrium, PyState, PyStateVec}; use feos_core::*; use feos_dft::adsorption::*; use feos_dft::interface::*; @@ -31,36 +28,18 @@ use feos_dft::*; use ndarray::{Array1, Array2, Array3, Array4}; use numpy::prelude::*; use numpy::{PyArray1, PyArray2, PyArray3, PyArray4}; -use pyo3::exceptions::{PyIndexError, PyValueError}; use pyo3::prelude::*; -#[cfg(feature = "estimator")] -use pyo3::wrap_pymodule; use quantity::*; -use std::collections::HashMap; use std::convert::TryInto; use std::sync::Arc; -use typenum::{Quot, P3}; - -type Functional = EquationOfState; +use typenum::Quot; #[pyclass(name = "HelmholtzEnergyFunctional")] #[derive(Clone)] -pub struct PyFunctionalVariant(pub Arc>); - -impl PyFunctionalVariant { - fn new(functional: DFT) -> Self - where - FunctionalVariant: From, - { - let functional: DFT = functional.into(); - let n = functional.components(); - let eos = functional.ideal_gas(IdealGasModel::NoModel(n)); - Self(Arc::new(eos)) - } -} +pub struct PyHelmholtzEnergyFunctional; #[pymethods] -impl PyFunctionalVariant { +impl PyHelmholtzEnergyFunctional { /// PC-SAFT Helmholtz energy functional. /// /// Parameters @@ -94,15 +73,20 @@ impl PyFunctionalVariant { max_iter_cross_assoc: usize, tol_cross_assoc: f64, dq_variant: DQVariants, - ) -> Self { + ) -> PyEquationOfState { + use super::eos::PyEquationOfState; + let options = PcSaftOptions { max_eta, max_iter_cross_assoc, tol_cross_assoc, dq_variant, }; - let func = PcSaftFunctional::with_options(parameters.0, fmt_version, options); - Self::new(func) + let func = Arc::new(ResidualModel::PcSaftFunctional( + PcSaftFunctional::with_options(parameters.0, fmt_version, options), + )); + let ideal_gas = Arc::new(IdealGasModel::NoModel(func.components())); + PyEquationOfState(Arc::new(EquationOfState::new(ideal_gas, func))) } /// (heterosegmented) group contribution PC-SAFT Helmholtz energy functional. @@ -135,14 +119,17 @@ impl PyFunctionalVariant { max_eta: f64, max_iter_cross_assoc: usize, tol_cross_assoc: f64, - ) -> Self { + ) -> PyEquationOfState { let options = GcPcSaftOptions { max_eta, max_iter_cross_assoc, tol_cross_assoc, }; - let func = GcPcSaftFunctional::with_options(parameters.0, fmt_version, options); - Self::new(func) + let func = Arc::new(ResidualModel::GcPcSaftFunctional( + GcPcSaftFunctional::with_options(parameters.0, fmt_version, options), + )); + let ideal_gas = Arc::new(IdealGasModel::NoModel(func.components())); + PyEquationOfState(Arc::new(EquationOfState::new(ideal_gas, func))) } /// PeTS Helmholtz energy functional without simplifications @@ -166,10 +153,19 @@ impl PyFunctionalVariant { signature = (parameters, fmt_version=FMTVersion::WhiteBear, max_eta=0.5), text_signature = "(parameters, fmt_version, max_eta=0.5)" )] - fn pets(parameters: PyPetsParameters, fmt_version: FMTVersion, max_eta: f64) -> Self { + fn pets( + parameters: PyPetsParameters, + fmt_version: FMTVersion, + max_eta: f64, + ) -> PyEquationOfState { let options = PetsOptions { max_eta }; - let func = PetsFunctional::with_options(parameters.0, fmt_version, options); - Self::new(func) + let func = Arc::new(ResidualModel::PetsFunctional(PetsFunctional::with_options( + parameters.0, + fmt_version, + options, + ))); + let ideal_gas = Arc::new(IdealGasModel::NoModel(func.components())); + PyEquationOfState(Arc::new(EquationOfState::new(ideal_gas, func))) } /// Helmholtz energy functional for hard sphere systems. @@ -185,9 +181,13 @@ impl PyFunctionalVariant { /// ------- /// HelmholtzEnergyFunctional #[staticmethod] - fn fmt(sigma: &Bound<'_, PyArray1>, fmt_version: FMTVersion) -> Self { - let func = FMTFunctional::new(&sigma.to_owned_array(), fmt_version); - Self::new(func) + fn fmt(sigma: &Bound<'_, PyArray1>, fmt_version: FMTVersion) -> PyEquationOfState { + let func = Arc::new(ResidualModel::FmtFunctional(FMTFunctional::new( + &sigma.to_owned_array(), + fmt_version, + ))); + let ideal_gas = Arc::new(IdealGasModel::NoModel(func.components())); + PyEquationOfState(Arc::new(EquationOfState::new(ideal_gas, func))) } /// SAFT-VRQ Mie Helmholtz energy functional. @@ -217,44 +217,32 @@ impl PyFunctionalVariant { fmt_version: FMTVersion, max_eta: f64, inc_nonadd_term: bool, - ) -> Self { + ) -> PyEquationOfState { let options = SaftVRQMieOptions { max_eta, inc_nonadd_term, }; - let func = SaftVRQMieFunctional::with_options(parameters.0, fmt_version, options); - Self::new(func) + let func = Arc::new(ResidualModel::SaftVRQMieFunctional( + SaftVRQMieFunctional::with_options(parameters.0, fmt_version, options), + )); + let ideal_gas = Arc::new(IdealGasModel::NoModel(func.components())); + PyEquationOfState(Arc::new(EquationOfState::new(ideal_gas, func))) } } -impl_equation_of_state!(PyFunctionalVariant); - -impl_state!(DFT, PyFunctionalVariant); -impl_phase_equilibrium!(DFT, PyFunctionalVariant); +impl_planar_interface!(EquationOfState); +impl_surface_tension_diagram!(EquationOfState); -impl_planar_interface!(Functional); -impl_surface_tension_diagram!(Functional); +impl_pore!(EquationOfState, PyEquationOfState); +impl_adsorption!(EquationOfState, PyEquationOfState); -impl_pore!(Functional, PyFunctionalVariant); -impl_adsorption!(Functional, PyFunctionalVariant); - -impl_pair_correlation!(Functional); -impl_solvation_profile!(Functional); - -#[cfg(feature = "estimator")] -impl_estimator!(DFT, PyFunctionalVariant); +impl_pair_correlation!(EquationOfState); +impl_solvation_profile!(EquationOfState); #[pymodule] pub fn dft(m: &Bound<'_, PyModule>) -> PyResult<()> { - m.add_class::()?; - m.add_class::()?; - - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; @@ -269,16 +257,5 @@ pub fn dft(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; m.add_class::()?; - #[cfg(feature = "estimator")] - m.add_wrapped(wrap_pymodule!(estimator_dft))?; - Ok(()) } - -#[cfg(feature = "estimator")] -#[pymodule] -pub fn estimator_dft(m: &Bound<'_, PyModule>) -> PyResult<()> { - m.add_class::()?; - m.add_class::()?; - m.add_class::() -} diff --git a/src/python/eos.rs b/src/python/eos.rs index c6411a42e..ea3e6cc77 100644 --- a/src/python/eos.rs +++ b/src/python/eos.rs @@ -1,4 +1,3 @@ -use crate::eos::ResidualModel; #[cfg(feature = "epcsaft")] use crate::epcsaft::python::PyElectrolytePcSaftParameters; #[cfg(feature = "epcsaft")] @@ -34,6 +33,7 @@ use crate::saftvrqmie::{SaftVRQMie, SaftVRQMieOptions}; use crate::uvtheory::python::PyUVTheoryParameters; #[cfg(feature = "uvtheory")] use crate::uvtheory::{Perturbation, UVTheory, UVTheoryOptions}; +use crate::ResidualModel; use super::dippr::PyDippr; use super::joback::PyJoback; diff --git a/src/python/mod.rs b/src/python/mod.rs index c60a2caf0..891383971 100644 --- a/src/python/mod.rs +++ b/src/python/mod.rs @@ -34,7 +34,6 @@ use dft::dft as dft_module; #[pymodule] pub fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add("__version__", env!("CARGO_PKG_VERSION"))?; - // m.add_wrapped(wrap_pymodule!(quantity_module))?; m.add_wrapped(wrap_pymodule!(eos_module))?; #[cfg(feature = "dft")] @@ -62,8 +61,6 @@ pub fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { set_path(m, "feos.eos.estimator", "eos.estimator_eos")?; #[cfg(feature = "dft")] set_path(m, "feos.dft", "dft")?; - #[cfg(all(feature = "dft", feature = "estimator"))] - set_path(m, "feos.dft.estimator", "dft.estimator_dft")?; set_path(m, "feos.joback", "joback")?; set_path(m, "feos.dippr", "dippr")?; set_path(m, "feos.cubic", "cubic")?; diff --git a/src/saftvrqmie/dft/mod.rs b/src/saftvrqmie/dft/mod.rs index dd81e02eb..24cfae237 100644 --- a/src/saftvrqmie/dft/mod.rs +++ b/src/saftvrqmie/dft/mod.rs @@ -3,14 +3,12 @@ use crate::saftvrqmie::eos::SaftVRQMieOptions; use crate::saftvrqmie::parameters::SaftVRQMieParameters; use dispersion::AttractiveFunctional; use feos_core::parameter::Parameter; -use feos_core::{Components, EosResult, Molarweight}; +use feos_core::{Components, EosResult, Molarweight, Residual, StateHD}; use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; -use feos_dft::{ - FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, -}; -use ndarray::{Array, Array1, Array2, ArrayView2, ScalarOperand}; +use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape}; +use ndarray::{Array, Array1, Array2, ScalarOperand}; use non_additive_hs::NonAddHardSphereFunctional; use num_dual::DualNum; use quantity::{MolarWeight, GRAM, MOL}; @@ -28,7 +26,7 @@ pub struct SaftVRQMieFunctional { } impl SaftVRQMieFunctional { - pub fn new(parameters: Arc) -> DFT { + pub fn new(parameters: Arc) -> Self { Self::with_options( parameters, FMTVersion::WhiteBear, @@ -36,7 +34,7 @@ impl SaftVRQMieFunctional { ) } - pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> DFT { + pub fn new_full(parameters: Arc, fmt_version: FMTVersion) -> Self { Self::with_options(parameters, fmt_version, SaftVRQMieOptions::default()) } @@ -44,12 +42,12 @@ impl SaftVRQMieFunctional { parameters: Arc, fmt_version: FMTVersion, saft_options: SaftVRQMieOptions, - ) -> DFT { - DFT(Self { + ) -> Self { + Self { parameters, fmt_version, options: saft_options, - }) + } } } @@ -64,19 +62,27 @@ impl Components for SaftVRQMieFunctional { self.fmt_version, self.options, ) - .0 } } -impl HelmholtzEnergyFunctional for SaftVRQMieFunctional { - type Contribution = SaftVRQMieFunctionalContribution; - +impl Residual for SaftVRQMieFunctional { fn compute_max_density(&self, moles: &Array1) -> f64 { self.options.max_eta * moles.sum() / (FRAC_PI_6 * &self.parameters.m * self.parameters.sigma.mapv(|v| v.powi(3)) * moles) .sum() } + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> Vec<(String, D)> { + self.evaluate_bulk(state) + } +} + +impl HelmholtzEnergyFunctional for SaftVRQMieFunctional { + type Contribution = SaftVRQMieFunctionalContribution; + fn contributions(&self) -> Box<(dyn Iterator)> { let mut contributions = Vec::with_capacity(3); diff --git a/tests/pcsaft/dft.rs b/tests/pcsaft/dft.rs index be1953899..772ada138 100644 --- a/tests/pcsaft/dft.rs +++ b/tests/pcsaft/dft.rs @@ -5,9 +5,9 @@ use feos::hard_sphere::FMTVersion; use feos::ideal_gas::Joback; use feos::pcsaft::{PcSaft, PcSaftFunctional, PcSaftParameters}; use feos_core::parameter::{IdentifierOption, Parameter}; -use feos_core::{Contributions, PhaseEquilibrium, State, Verbosity}; +use feos_core::{Contributions, EquationOfState, PhaseEquilibrium, State, Verbosity}; use feos_dft::interface::PlanarInterface; -use feos_dft::DFTSolver; +use feos_dft::{DFTSolver, PdgtFunctionalProperties}; use ndarray::{arr1, Axis}; use quantity::*; use std::error::Error; @@ -337,7 +337,10 @@ fn test_entropy_bulk_values() -> Result<(), Box> { None, IdentifierOption::Name, )?; - let func = Arc::new(PcSaftFunctional::new(Arc::new(params)).ideal_gas(joback)); + let func = Arc::new(EquationOfState::new( + Arc::new(joback), + Arc::new(PcSaftFunctional::new(Arc::new(params))), + )); let vle = PhaseEquilibrium::pure(&func, 350.0 * KELVIN, None, Default::default())?; let profile = PlanarInterface::from_pdgt(&vle, 2048, false)?.solve(None)?; let s_res = profile.profile.entropy_density(Contributions::Residual)?; From 9f7b33611b4e6d7238b8883a36dde61a34bfa323 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 13 Mar 2025 11:57:42 +0100 Subject: [PATCH 2/3] minor fix --- feos-derive/src/lib.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/feos-derive/src/lib.rs b/feos-derive/src/lib.rs index 2304446b4..f9a374d7e 100644 --- a/feos-derive/src/lib.rs +++ b/feos-derive/src/lib.rs @@ -84,7 +84,7 @@ pub fn derive_residual(input: TokenStream) -> TokenStream { .into() } -#[proc_macro_derive(HelmholtzEnergyFunctional, attributes(implement_dft))] +#[proc_macro_derive(HelmholtzEnergyFunctional, attributes(implement))] pub fn derive_helmholtz_energy_functional(input: TokenStream) -> TokenStream { let input = parse_macro_input!(input as DeriveInput); expand_helmholtz_energy_functional(input) From 97d254c972eaa1acf43351b4bf1a25c2895b01eb Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 13 Mar 2025 12:01:53 +0100 Subject: [PATCH 3/3] fix conditional compiling --- src/eos.rs | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/eos.rs b/src/eos.rs index e13ad5a17..3dbf2ebe7 100644 --- a/src/eos.rs +++ b/src/eos.rs @@ -54,22 +54,23 @@ pub enum ResidualModel { UVTheory(crate::uvtheory::UVTheory), // Helmholtz energy functionals - #[cfg(feature = "dft")] + #[cfg(all(feature = "dft", feature = "pcsaft"))] #[implement(molar_weight, functional, fluid_parameters, pair_potential)] PcSaftFunctional(crate::pcsaft::PcSaftFunctional), - #[cfg(feature = "gc_pcsaft")] + #[cfg(all(feature = "dft", feature = "gc_pcsaft"))] #[implement(molar_weight, functional, fluid_parameters, bond_lengths)] GcPcSaftFunctional(crate::gc_pcsaft::GcPcSaftFunctional), - #[cfg(feature = "pets")] + #[cfg(all(feature = "dft", feature = "pets"))] #[implement(molar_weight, functional, fluid_parameters, pair_potential)] PetsFunctional(crate::pets::PetsFunctional), + #[cfg(feature = "dft")] #[implement(functional, fluid_parameters, pair_potential)] FmtFunctional(crate::hard_sphere::FMTFunctional), - #[cfg(feature = "saftvrqmie")] + #[cfg(all(feature = "dft", feature = "saftvrqmie"))] #[implement(molar_weight, functional, fluid_parameters, pair_potential)] SaftVRQMieFunctional(crate::saftvrqmie::SaftVRQMieFunctional), }