diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs index c0ea3d357..1b178ef51 100644 --- a/feos-core/src/cubic.rs +++ b/feos-core/src/cubic.rs @@ -8,7 +8,7 @@ use crate::equation_of_state::{Components, Residual}; use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord}; use crate::si::{MolarWeight, GRAM, MOL}; use crate::state::StateHD; -use ndarray::{Array1, Array2}; +use ndarray::{Array1, Array2, ScalarOperand}; use num_dual::DualNum; use serde::{Deserialize, Serialize}; use std::f64::consts::SQRT_2; @@ -146,16 +146,6 @@ impl Parameter for PengRobinsonParameters { } } -struct PengRobinsonContribution { - parameters: Arc, -} - -impl fmt::Display for PengRobinsonContribution { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Peng Robinson") - } -} - /// A simple version of the Peng-Robinson equation of state. pub struct PengRobinson { /// Parameters @@ -215,7 +205,7 @@ impl Residual for PengRobinson { * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln()) } - fn residual_helmholtz_energy_contributions + Copy>( + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( &self, state: &StateHD, ) -> Vec<(String, D)> { diff --git a/feos-core/src/equation_of_state/ideal_gas.rs b/feos-core/src/equation_of_state/ideal_gas.rs index 99939d034..116c8ee11 100644 --- a/feos-core/src/equation_of_state/ideal_gas.rs +++ b/feos-core/src/equation_of_state/ideal_gas.rs @@ -2,12 +2,13 @@ use super::Components; use crate::StateHD; use ndarray::Array1; use num_dual::DualNum; -use std::fmt::Display; /// Ideal gas Helmholtz energy contribution. -pub trait IdealGas: Components + Sync + Send + Display { +pub trait IdealGas: Components + Sync + Send { fn ln_lambda3 + Copy>(&self, temperature: D) -> Array1; + fn ideal_gas_model(&self) -> String; + /// Evaluate the ideal gas Helmholtz energy contribution for a given state. /// /// In some cases it could be advantageous to overwrite this diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs index 54afdee50..139031e67 100644 --- a/feos-core/src/equation_of_state/mod.rs +++ b/feos-core/src/equation_of_state/mod.rs @@ -2,8 +2,8 @@ use crate::{ si::{Diffusivity, MolarWeight, Moles, Temperature, ThermalConductivity, Viscosity, Volume}, EosResult, }; -use ndarray::Array1; -use std::{fmt::Display, sync::Arc}; +use ndarray::{Array1, ScalarOperand}; +use std::sync::Arc; mod ideal_gas; mod residual; @@ -52,12 +52,6 @@ impl EquationOfState { } } -impl Display for EquationOfState { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "{}", self.ideal_gas) - } -} - impl Components for EquationOfState { fn components(&self) -> usize { assert_eq!( @@ -80,6 +74,10 @@ impl IdealGas for EquationOfState + Copy>(&self, temperature: D) -> Array1 { self.ideal_gas.ln_lambda3(temperature) } + + fn ideal_gas_model(&self) -> String { + self.ideal_gas.ideal_gas_model() + } } impl Residual for EquationOfState { @@ -87,14 +85,7 @@ impl Residual for EquationOfState { self.residual.compute_max_density(moles) } - fn residual_helmholtz_energy + Copy>( - &self, - state: &crate::StateHD, - ) -> D { - self.residual.residual_helmholtz_energy(state) - } - - fn residual_helmholtz_energy_contributions + Copy>( + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( &self, state: &crate::StateHD, ) -> Vec<(String, D)> { diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs index 84309a7d4..2ea004b63 100644 --- a/feos-core/src/equation_of_state/residual.rs +++ b/feos-core/src/equation_of_state/residual.rs @@ -3,6 +3,7 @@ use crate::si::*; use crate::StateHD; use crate::{EosError, EosResult}; use ndarray::prelude::*; +use ndarray::ScalarOperand; use num_dual::*; use num_traits::{One, Zero}; use std::ops::Div; @@ -24,13 +25,16 @@ pub trait Residual: Components + Send + Sync { /// Evaluate the reduced Helmholtz energy of each individual contribution /// and return them together with a string representation of the contribution. - fn residual_helmholtz_energy_contributions + Copy>( + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( &self, state: &StateHD, ) -> Vec<(String, D)>; /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$. - fn residual_helmholtz_energy + Copy>(&self, state: &StateHD) -> D { + fn residual_helmholtz_energy + Copy + ScalarOperand>( + &self, + state: &StateHD, + ) -> D { self.residual_helmholtz_energy_contributions(state) .iter() .fold(D::zero(), |acc, (_, a)| acc + a) diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs index ae5c7f0f7..7c850261b 100644 --- a/feos-core/src/lib.rs +++ b/feos-core/src/lib.rs @@ -126,18 +126,11 @@ mod tests { use crate::EosResult; use crate::StateBuilder; use approx::*; - use std::fmt::Display; use std::sync::Arc; // Only to be able to instantiate an `EquationOfState` struct NoIdealGas; - impl Display for NoIdealGas { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "NoIdealGas") - } - } - impl Components for NoIdealGas { fn components(&self) -> usize { 1 @@ -155,6 +148,10 @@ mod tests { ) -> ndarray::prelude::Array1 { unreachable!() } + + fn ideal_gas_model(&self) -> String { + "NoIdealGas".into() + } } fn pure_record_vec() -> Vec> { diff --git a/feos-core/src/parameter/mod.rs b/feos-core/src/parameter/mod.rs index 6daca05cf..8633e7825 100644 --- a/feos-core/src/parameter/mod.rs +++ b/feos-core/src/parameter/mod.rs @@ -84,7 +84,7 @@ where /// `pure_records`, the `Default` implementation of Self::Binary is used. #[allow(clippy::expect_fun_call)] fn binary_matrix_from_records( - pure_records: &Vec>, + pure_records: &[PureRecord], binary_records: &[BinaryRecord], identifier_option: IdentifierOption, ) -> Option> { diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs index 4990f84d7..7973bbcbd 100644 --- a/feos-core/src/state/properties.rs +++ b/feos-core/src/state/properties.rs @@ -257,7 +257,7 @@ impl State { let contributions = self.eos.residual_helmholtz_energy_contributions(&new_state); let mut res = Vec::with_capacity(contributions.len() + 1); res.push(( - self.eos.to_string(), + self.eos.ideal_gas_model(), MolarEnergy::from_reduced( (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature).eps, ), diff --git a/feos-core/tests/parameters.rs b/feos-core/tests/parameters.rs index d1fe0bfd8..1498b7028 100644 --- a/feos-core/tests/parameters.rs +++ b/feos-core/tests/parameters.rs @@ -82,7 +82,7 @@ fn from_records() -> Result<(), ParameterError> { } ] "#; - let pure_records = serde_json::from_str(pr_json).expect("Unable to parse json."); + let pure_records: Vec<_> = serde_json::from_str(pr_json).expect("Unable to parse json."); let binary_records: Vec<_> = serde_json::from_str(br_json).expect("Unable to parse json."); let binary_matrix = MyParameter::binary_matrix_from_records( &pure_records, @@ -201,7 +201,7 @@ fn from_records_missing_binary() -> Result<(), ParameterError> { } ] "#; - let pure_records = serde_json::from_str(pr_json).expect("Unable to parse json."); + let pure_records: Vec<_> = serde_json::from_str(pr_json).expect("Unable to parse json."); let binary_records: Vec<_> = serde_json::from_str(br_json).expect("Unable to parse json."); let binary_matrix = MyParameter::binary_matrix_from_records( &pure_records, @@ -266,7 +266,7 @@ fn from_records_correct_binary_order() -> Result<(), ParameterError> { } ] "#; - let pure_records = serde_json::from_str(pr_json).expect("Unable to parse json."); + let pure_records: Vec<_> = serde_json::from_str(pr_json).expect("Unable to parse json."); let binary_records: Vec<_> = serde_json::from_str(br_json).expect("Unable to parse json."); let binary_matrix = MyParameter::binary_matrix_from_records( &pure_records, diff --git a/feos-derive/src/dft.rs b/feos-derive/src/dft.rs index 5eedf4ee5..0b5772e99 100644 --- a/feos-derive/src/dft.rs +++ b/feos-derive/src/dft.rs @@ -95,7 +95,7 @@ fn impl_helmholtz_energy_functional( let contributions = variants.iter().map(|v| { let name = &v.ident; quote! { - Self::#name(functional) => functional.contributions() + Self::#name(functional) => Box::new(functional.contributions().map(FunctionalContributionVariant::from)) } }); @@ -121,6 +121,7 @@ fn impl_helmholtz_energy_functional( Ok(quote! { impl HelmholtzEnergyFunctional for FunctionalVariant { + type Contribution = FunctionalContributionVariant; fn molecule_shape(&self) -> MoleculeShape { match self { #(#molecule_shape,)* @@ -131,7 +132,7 @@ fn impl_helmholtz_energy_functional( #(#compute_max_density,)* } } - fn contributions(&self) -> &[Box] { + fn contributions(&self) -> Box> { match self { #(#contributions,)* } diff --git a/feos-derive/src/functional_contribution.rs b/feos-derive/src/functional_contribution.rs new file mode 100644 index 000000000..0f6a25487 --- /dev/null +++ b/feos-derive/src/functional_contribution.rs @@ -0,0 +1,115 @@ +use quote::quote; +use syn::{DeriveInput, Ident}; + +pub(crate) fn expand_functional_contribution( + input: DeriveInput, +) -> syn::Result { + let ident = input.ident; + let variants = match input.data { + syn::Data::Enum(syn::DataEnum { ref variants, .. }) => variants, + _ => panic!("this derive macro only works on enums"), + }; + + let functional_contribution = impl_functional_contribution(&ident, variants); + let display = impl_display(&ident, variants); + let from = impl_from(&ident, variants); + Ok(quote! { + #functional_contribution + #display + #from + }) +} + +fn impl_functional_contribution( + ident: &Ident, + variants: &syn::punctuated::Punctuated, +) -> proc_macro2::TokenStream { + let weight_functions = variants.iter().map(|v| { + let name = &v.ident; + quote! { + Self::#name(functional_contribution) => functional_contribution.weight_functions(temperature) + } + }); + let weight_functions_pdgt = variants.iter().map(|v| { + let name = &v.ident; + quote! { + Self::#name(functional_contribution) => functional_contribution.weight_functions_pdgt(temperature) + } + }); + let helmholtz_energy_density = variants.iter().map(|v| { + let name = &v.ident; + quote! { + Self::#name(functional_contribution) => functional_contribution.helmholtz_energy_density(temperature, weighted_densities) + } + }); + + quote! { + impl FunctionalContribution for #ident { + fn weight_functions + Copy+ScalarOperand>(&self, temperature: N) -> WeightFunctionInfo { + match self { + #(#weight_functions,)* + } + } + fn weight_functions_pdgt + Copy+ScalarOperand>(&self, temperature: N) -> WeightFunctionInfo { + match self { + #(#weight_functions_pdgt,)* + } + } + fn helmholtz_energy_density + Copy+ScalarOperand>( + &self, + temperature: N, + weighted_densities: ArrayView2, + ) -> EosResult> { + match self { + #(#helmholtz_energy_density,)* + } + } + } + } +} + +fn impl_display( + ident: &Ident, + variants: &syn::punctuated::Punctuated, +) -> proc_macro2::TokenStream { + let fmt = variants.iter().map(|v| { + let name = &v.ident; + quote! { + Self::#name(functional_contribution) => functional_contribution.fmt(f) + } + }); + + quote! { + impl std::fmt::Display for #ident { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + #(#fmt,)* + } + } + } + } +} + +fn impl_from( + ident: &Ident, + variants: &syn::punctuated::Punctuated, +) -> proc_macro2::TokenStream { + let from = variants.iter().map(|v| { + let name = &v.ident; + let syn::Fields::Unnamed(syn::FieldsUnnamed { unnamed, .. }) = &v.fields else { + panic!("All variants must be tuple structs!") + }; + let inner = &unnamed.first().unwrap().ty; + quote! { + impl From<#inner> for #ident { + fn from(variant: #inner) -> Self { + Self::#name(variant) + } + } + } + }); + + quote! { + #(#from)* + } +} diff --git a/feos-derive/src/ideal_gas.rs b/feos-derive/src/ideal_gas.rs index 86a3bbc06..e03cf959d 100644 --- a/feos-derive/src/ideal_gas.rs +++ b/feos-derive/src/ideal_gas.rs @@ -8,10 +8,8 @@ pub(crate) fn expand_ideal_gas(input: DeriveInput) -> syn::Result panic!("No ideal gas model initialized!") + } + } else { + quote! { + Self::#name(ideal_gas) => ideal_gas.ideal_gas_model() + } + } + }); quote! { impl IdealGas for IdealGasModel { fn ln_lambda3 + Copy>(&self, temperature: D) -> Array1 { @@ -37,26 +47,11 @@ fn impl_ideal_gas( #(#ln_lambda3,)* } } - } - } -} -fn impl_display( - variants: &syn::punctuated::Punctuated, -) -> proc_macro2::TokenStream { - let string = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(ideal_gas) => ideal_gas.to_string() - } - }); - quote! { - impl fmt::Display for IdealGasModel { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let s = match self { + fn ideal_gas_model(&self) -> String { + match self { #(#string,)* - }; - write!(f, "{}", s) + } } } } diff --git a/feos-derive/src/lib.rs b/feos-derive/src/lib.rs index e91a82b76..06ae18443 100644 --- a/feos-derive/src/lib.rs +++ b/feos-derive/src/lib.rs @@ -3,6 +3,7 @@ //! the boilerplate for the EquationOfState and HelmholtzEnergyFunctional traits. use components::expand_components; use dft::expand_helmholtz_energy_functional; +use functional_contribution::expand_functional_contribution; use ideal_gas::expand_ideal_gas; use proc_macro::TokenStream; use residual::expand_residual; @@ -10,6 +11,7 @@ use syn::{parse_macro_input, DeriveInput}; mod components; mod dft; +mod functional_contribution; mod ideal_gas; mod residual; @@ -78,3 +80,11 @@ pub fn derive_helmholtz_energy_functional(input: TokenStream) -> TokenStream { .unwrap_or_else(syn::Error::into_compile_error) .into() } + +#[proc_macro_derive(FunctionalContribution, attributes(implement))] +pub fn derive_functional_contribution(input: TokenStream) -> TokenStream { + let input = parse_macro_input!(input as DeriveInput); + expand_functional_contribution(input) + .unwrap_or_else(syn::Error::into_compile_error) + .into() +} diff --git a/feos-derive/src/residual.rs b/feos-derive/src/residual.rs index b68e4da43..4604e87a1 100644 --- a/feos-derive/src/residual.rs +++ b/feos-derive/src/residual.rs @@ -28,12 +28,6 @@ fn impl_residual( Self::#name(residual) => residual.compute_max_density(moles) } }); - let residual_helmholtz_energy = variants.iter().map(|v| { - let name = &v.ident; - quote! { - Self::#name(residual) => residual.residual_helmholtz_energy(state) - } - }); let residual_helmholtz_energy_contributions = variants.iter().map(|v| { let name = &v.ident; quote! { @@ -54,12 +48,7 @@ fn impl_residual( #(#compute_max_density,)* } } - fn residual_helmholtz_energy + Copy>(&self, state: &StateHD) -> D { - match self { - #(#residual_helmholtz_energy,)* - } - } - fn residual_helmholtz_energy_contributions + Copy>(&self, state: &StateHD) -> Vec<(String, D)> { + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>(&self, state: &StateHD) -> Vec<(String, D)> { match self { #(#residual_helmholtz_energy_contributions,)* } diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index bfb0a4fba..7efe91ec6 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -5,6 +5,7 @@ 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::si::{ Density, Dimensionless, Energy, Length, MolarEnergy, MolarWeight, Temperature, Volume, KELVIN, }; @@ -13,6 +14,8 @@ use ndarray::prelude::*; use ndarray::Axis as Axis_nd; use ndarray::RemoveAxis; use num_dual::linalg::LU; +use num_dual::DualNum; +use std::fmt::Display; use std::sync::Arc; const POTENTIAL_OFFSET: f64 = 2.0; @@ -283,8 +286,10 @@ impl Components for Helium { } impl HelmholtzEnergyFunctional for Helium { - fn contributions(&self) -> &[Box] { - &[] + type Contribution = HeliumContribution; + + fn contributions(&self) -> Box> { + Box::new([].into_iter()) } fn compute_max_density(&self, _: &Array1) -> f64 { @@ -309,3 +314,25 @@ impl FluidParameters for Helium { &self.sigma } } + +struct HeliumContribution; + +impl FunctionalContribution for HeliumContribution { + fn weight_functions + Copy>(&self, _: N) -> WeightFunctionInfo { + unreachable!() + } + + fn helmholtz_energy_density + Copy>( + &self, + _: N, + _: ArrayView2, + ) -> EosResult> { + unreachable!() + } +} + +impl Display for HeliumContribution { + fn fmt(&self, _: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + unreachable!() + } +} diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 913406f10..e825b12bd 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -5,10 +5,7 @@ use crate::ideal_chain_contribution::IdealChainContribution; use crate::solvation::PairPotential; use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use feos_core::si::MolarWeight; -use feos_core::{ - Components, DeBroglieWavelength, EosResult, EquationOfState, HelmholtzEnergy, - HelmholtzEnergyDual, IdealGas, Residual, StateHD, -}; +use feos_core::{Components, EosResult, EquationOfState, IdealGas, Residual, StateHD}; use ndarray::*; use num_dual::*; use petgraph::graph::{Graph, UnGraph}; @@ -21,7 +18,9 @@ use std::sync::Arc; impl HelmholtzEnergyFunctional for EquationOfState { - fn contributions(&self) -> &[Box] { + type Contribution = F::Contribution; + + fn contributions(&self) -> Box> { self.residual.contributions() } @@ -99,43 +98,18 @@ impl Residual for DFT { self.0.compute_max_density(moles) } - fn contributions(&self) -> &[Box] { - unreachable!() - } - fn molar_weight(&self) -> MolarWeight> { self.0.molar_weight() } - fn evaluate_residual + Copy>(&self, state: &StateHD) -> D - where - dyn HelmholtzEnergy: HelmholtzEnergyDual, - { - self.0 - .contributions() - .iter() - .map(|c| (c as &dyn HelmholtzEnergy).helmholtz_energy(state)) - .sum::() - + self.ideal_chain_contribution().helmholtz_energy(state) - } - - fn evaluate_residual_contributions + Copy>( + fn residual_helmholtz_energy_contributions + Copy + ScalarOperand>( &self, state: &StateHD, - ) -> Vec<(String, D)> - where - dyn HelmholtzEnergy: HelmholtzEnergyDual, - { + ) -> Vec<(String, D)> { let mut res: Vec<(String, D)> = self .0 .contributions() - .iter() - .map(|c| { - ( - c.to_string(), - (c as &dyn HelmholtzEnergy).helmholtz_energy(state), - ) - }) + .map(|c| (c.to_string(), c.helmholtz_energy(state))) .collect(); res.push(( self.ideal_chain_contribution().to_string(), @@ -146,7 +120,11 @@ impl Residual for DFT { } impl IdealGas for DFT { - fn ideal_gas_model(&self) -> &dyn DeBroglieWavelength { + fn ln_lambda3 + Copy>(&self, temperature: D) -> Array1 { + self.0.ln_lambda3(temperature) + } + + fn ideal_gas_model(&self) -> String { self.0.ideal_gas_model() } } @@ -165,8 +143,10 @@ pub enum MoleculeShape<'a> { /// A general Helmholtz energy functional. pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { + type Contribution: FunctionalContribution; + /// Return a slice of [FunctionalContribution]s. - fn contributions(&self) -> &[Box]; + fn contributions(&self) -> Box>; /// Return the shape of the molecules and the necessary specifications. fn molecule_shape(&self) -> MoleculeShape; @@ -191,7 +171,6 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { fn weight_functions(&self, temperature: f64) -> Vec> { self.contributions() - .iter() .map(|c| c.weight_functions(temperature)) .collect() } @@ -232,9 +211,9 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { { let weighted_densities = convolver.weighted_densities(density); let contributions = self.contributions(); - let mut partial_derivatives = Vec::with_capacity(contributions.len()); + let mut partial_derivatives = Vec::new(); let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0))); - for (c, wd) in contributions.iter().zip(weighted_densities) { + for (c, wd) in contributions.zip(weighted_densities) { let nwd = wd.shape()[0]; let ngrid = wd.len() / nwd; let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); @@ -269,9 +248,9 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { let density_dual = density.mapv(Dual64::from); let weighted_densities = convolver.weighted_densities(&density_dual); let contributions = self.contributions(); - let mut partial_derivatives = Vec::with_capacity(contributions.len()); + let mut partial_derivatives = Vec::new(); let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0))); - for (c, wd) in contributions.iter().zip(weighted_densities) { + for (c, wd) in contributions.zip(weighted_densities) { let nwd = wd.shape()[0]; let ngrid = wd.len() / nwd; let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index 89bbf8499..b714fa48f 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -1,97 +1,51 @@ use crate::weight_functions::WeightFunctionInfo; -use feos_core::{EosResult, HelmholtzEnergyDual, StateHD}; +use feos_core::{EosResult, StateHD}; use ndarray::prelude::*; -use ndarray::RemoveAxis; +use ndarray::{RemoveAxis, ScalarOperand}; use num_dual::*; use num_traits::{One, Zero}; use std::fmt::Display; -macro_rules! impl_helmholtz_energy { - ($number:ty) => { - impl HelmholtzEnergyDual<$number> for Box { - fn helmholtz_energy(&self, state: &StateHD<$number>) -> $number { - // calculate weight functions - let weight_functions = self.weight_functions(state.temperature); - - // calculate segment density - let density = weight_functions - .component_index - .mapv(|c| state.partial_density[c]); - - // calculate weighted density and Helmholtz energy - let weight_constants = weight_functions.weight_constants(Zero::zero(), 0); - let weighted_densities = weight_constants.dot(&density).insert_axis(Axis(1)); - self.calculate_helmholtz_energy_density( - state.temperature, - weighted_densities.view(), - ) - .unwrap()[0] - * state.volume - } - } - }; -} - -impl_helmholtz_energy!(f64); -impl_helmholtz_energy!(Dual64); -impl_helmholtz_energy!(Dual, f64>); -impl_helmholtz_energy!(HyperDual64); -impl_helmholtz_energy!(Dual2_64); -impl_helmholtz_energy!(Dual3_64); -impl_helmholtz_energy!(HyperDual); -impl_helmholtz_energy!(HyperDual, f64>); -impl_helmholtz_energy!(HyperDual, f64>); -impl_helmholtz_energy!(Dual2); -impl_helmholtz_energy!(Dual3); -impl_helmholtz_energy!(Dual3, f64>); -impl_helmholtz_energy!(Dual3, f64>); - -/// Individual functional contribution that can -/// be evaluated using generalized (hyper) dual numbers. -/// -/// This trait needs to be implemented generically or for -/// the specific types in the supertraits of [FunctionalContribution] -/// so that the implementor can be used as a functional -/// contribution in the Helmholtz energy functional. -pub trait FunctionalContributionDual>: Display { +/// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers. +pub trait FunctionalContribution: Display + Sync + Send { /// Return the weight functions required in this contribution. - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo; + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo; + /// Overwrite this if the weight functions in pDGT are different than for DFT. - fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { self.weight_functions(temperature) } /// Return the Helmholtz energy density for the given temperature and weighted densities. - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, ) -> EosResult>; -} -/// Object safe version of the [FunctionalContributionDual] trait. -/// -/// The trait is implemented automatically for every struct that implements -/// the supertraits. -pub trait FunctionalContribution: - FunctionalContributionDual - + FunctionalContributionDual - + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual - + FunctionalContributionDual - + FunctionalContributionDual - + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual> - + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> - + Display - + Sync - + Send -{ + fn helmholtz_energy + Copy + ScalarOperand>(&self, state: &StateHD) -> N { + // calculate weight functions + let weight_functions = self.weight_functions(state.temperature); + + // calculate segment density + let density = weight_functions + .component_index + .mapv(|c| state.partial_density[c]); + + // calculate weighted density and Helmholtz energy + let weight_constants = weight_functions.weight_constants(Zero::zero(), 0); + let weighted_densities = weight_constants.dot(&density).insert_axis(Axis(1)); + self.helmholtz_energy_density(state.temperature, weighted_densities.view()) + .unwrap()[0] + * state.volume + } + fn first_partial_derivatives( &self, temperature: f64, @@ -105,7 +59,7 @@ pub trait FunctionalContribution: for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps = 1.0); - phi = self.calculate_helmholtz_energy_density(t, wd.view())?; + phi = self.helmholtz_energy_density(t, wd.view())?; first_partial_derivative .index_axis_mut(Axis(0), i) .assign(&phi.mapv(|p| p.eps)); @@ -129,7 +83,7 @@ pub trait FunctionalContribution: for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i) .map_inplace(|x| x.eps = Dual::one()); - phi = self.calculate_helmholtz_energy_density(t, wd.view())?; + phi = self.helmholtz_energy_density(t, wd.view())?; first_partial_derivative .index_axis_mut(Axis(0), i) .assign(&phi.mapv(|p| p.eps)); @@ -156,7 +110,7 @@ pub trait FunctionalContribution: wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps1 = 1.0); for j in 0..=i { wd.index_axis_mut(Axis(0), j).map_inplace(|x| x.eps2 = 1.0); - phi = self.calculate_helmholtz_energy_density(t, wd.view())?; + phi = self.helmholtz_energy_density(t, wd.view())?; let p = phi.mapv(|p| p.eps1eps2); second_partial_derivative .index_axis_mut(Axis(0), i) @@ -179,24 +133,3 @@ pub trait FunctionalContribution: Ok(()) } } - -impl FunctionalContribution for T where - T: FunctionalContributionDual - + FunctionalContributionDual - + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual - + FunctionalContributionDual - + FunctionalContributionDual - + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual> - + FunctionalContributionDual> - + FunctionalContributionDual, f64>> - + FunctionalContributionDual, f64>> - + Display - + Sync - + Send -{ -} diff --git a/feos-dft/src/ideal_chain_contribution.rs b/feos-dft/src/ideal_chain_contribution.rs index d42e68a00..5965b2553 100644 --- a/feos-dft/src/ideal_chain_contribution.rs +++ b/feos-dft/src/ideal_chain_contribution.rs @@ -1,5 +1,5 @@ use feos_core::si::{Density, Pressure, Temperature}; -use feos_core::{EosResult, HelmholtzEnergyDual, StateHD}; +use feos_core::{EosResult, StateHD}; use ndarray::*; use num_dual::DualNum; use std::fmt; @@ -19,8 +19,14 @@ impl IdealChainContribution { } } -impl + Copy> HelmholtzEnergyDual for IdealChainContribution { - fn helmholtz_energy(&self, state: &StateHD) -> D { +impl fmt::Display for IdealChainContribution { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "Ideal chain") + } +} + +impl IdealChainContribution { + pub fn helmholtz_energy + Copy>(&self, state: &StateHD) -> D { let segments = self.component_index.len(); if self.component_index[segments - 1] + 1 != segments { return D::zero(); @@ -36,16 +42,8 @@ impl + Copy> HelmholtzEnergyDual for IdealChainContribution { .sum() * state.volume } -} -impl fmt::Display for IdealChainContribution { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Ideal chain") - } -} - -impl IdealChainContribution { - pub fn calculate_helmholtz_energy_density( + pub fn helmholtz_energy_density( &self, density: &Array, ) -> EosResult> @@ -61,7 +59,7 @@ impl IdealChainContribution { Ok(phi) } - pub fn helmholtz_energy_density( + pub fn helmholtz_energy_density_units( &self, temperature: Temperature, density: &Density>, @@ -73,7 +71,7 @@ impl IdealChainContribution { let rho = density.to_reduced(); let t = temperature.to_reduced(); Ok(Pressure::from_reduced( - self.calculate_helmholtz_energy_density(&rho)? * t, + self.helmholtz_energy_density(&rho)? * t, )) } } diff --git a/feos-dft/src/lib.rs b/feos-dft/src/lib.rs index 2962fbeea..58f65a5a3 100644 --- a/feos-dft/src/lib.rs +++ b/feos-dft/src/lib.rs @@ -18,7 +18,7 @@ mod weight_functions; pub use convolver::{Convolver, ConvolverFFT}; pub use functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; -pub use functional_contribution::{FunctionalContribution, FunctionalContributionDual}; +pub use functional_contribution::FunctionalContribution; pub use geometry::{Axis, Geometry, Grid}; pub use profile::{DFTProfile, DFTSpecification, DFTSpecifications}; pub use solver::{DFTSolver, DFTSolverLog}; diff --git a/feos-dft/src/pdgt.rs b/feos-dft/src/pdgt.rs index 4cd836bef..162f58440 100644 --- a/feos-dft/src/pdgt.rs +++ b/feos-dft/src/pdgt.rs @@ -12,18 +12,18 @@ use std::ops::{Add, AddAssign, Sub}; use typenum::{Diff, Sum, P2}; type _InfluenceParameter = Diff, _Density>; -type InfluenceParameter = Quantity; +pub type InfluenceParameter = Quantity; impl WeightFunctionInfo { - fn pdgt_weight_constants(&self) -> (Array2, Array2, Array2) { + pub fn pdgt_weight_constants(&self) -> (Array2, Array2, Array2) { let k = Dual2_64::from(0.0).derivative(); let w = self.weight_constants(k, 1); (w.mapv(|w| w.re), w.mapv(|w| -w.v1), w.mapv(|w| -0.5 * w.v2)) } } -impl dyn FunctionalContribution { - pub fn pdgt_properties( +trait PdgtProperties: FunctionalContribution { + fn pdgt_properties( &self, temperature: f64, density: &Array2, @@ -116,7 +116,7 @@ impl dyn FunctionalContribution { } #[allow(clippy::type_complexity)] - pub fn influence_diagonal( + fn influence_diagonal( &self, temperature: Temperature, density: &Density>, @@ -141,6 +141,8 @@ impl dyn FunctionalContribution { } } +impl PdgtProperties for T {} + impl DFT { pub fn solve_pdgt( &self, @@ -172,7 +174,7 @@ impl DFT { } delta_omega += &self .ideal_chain_contribution() - .helmholtz_energy_density::(vle.vapor().temperature, &density)?; + .helmholtz_energy_density_units::(vle.vapor().temperature, &density)?; // calculate excess grand potential density let mu_res = vle.vapor().residual_chemical_potential(); diff --git a/feos-dft/src/profile/properties.rs b/feos-dft/src/profile/properties.rs index 9afcdb9aa..503bb1481 100644 --- a/feos-dft/src/profile/properties.rs +++ b/feos-dft/src/profile/properties.rs @@ -1,7 +1,7 @@ #![allow(type_alias_bounds)] use super::DFTProfile; use crate::convolver::{BulkConvolver, Convolver}; -use crate::functional_contribution::{FunctionalContribution, FunctionalContributionDual}; +use crate::functional_contribution::FunctionalContribution; use crate::{ConvolverFFT, DFTSolverLog, HelmholtzEnergyFunctional, WeightFunctionInfo}; use feos_core::si::{ Density, Energy, Entropy, EntropyDensity, MolarEnergy, Moles, Pressure, Quantity, Temperature, @@ -79,7 +79,6 @@ where ) -> EosResult> where N: DualNum + Copy + ScalarOperand, - dyn FunctionalContribution: FunctionalContributionDual, { let density_dual = density.mapv(N::from); let weighted_densities = convolver.weighted_densities(&density_dual); @@ -87,15 +86,15 @@ where let mut helmholtz_energy_density: Array = self .dft .ideal_chain_contribution() - .calculate_helmholtz_energy_density(&density.mapv(N::from))?; - for (c, wd) in functional_contributions.iter().zip(weighted_densities) { + .helmholtz_energy_density(&density.mapv(N::from))?; + for (c, wd) in functional_contributions.zip(weighted_densities) { let nwd = wd.shape()[0]; let ngrid = wd.len() / nwd; helmholtz_energy_density .view_mut() .into_shape(ngrid) .unwrap() - .add_assign(&c.calculate_helmholtz_energy_density( + .add_assign(&c.helmholtz_energy_density( temperature, wd.into_shape((nwd, ngrid)).unwrap().view(), )?); @@ -112,7 +111,6 @@ where let temperature_dual = Dual64::from(temperature).derivative(); let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions - .iter() .map(|c| c.weight_functions(temperature_dual)) .collect(); let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None); @@ -139,19 +137,18 @@ where let temperature_dual = Dual64::from(temperature).derivative(); let weighted_densities = convolver.weighted_densities(&density_dual); let functional_contributions = self.dft.contributions(); - let mut helmholtz_energy_density: Vec> = - Vec::with_capacity(functional_contributions.len() + 1); + let mut helmholtz_energy_density: Vec> = Vec::new(); helmholtz_energy_density.push( self.dft .ideal_chain_contribution() - .calculate_helmholtz_energy_density(&density.mapv(Dual64::from))?, + .helmholtz_energy_density(&density.mapv(Dual64::from))?, ); - for (c, wd) in functional_contributions.iter().zip(weighted_densities) { + for (c, wd) in functional_contributions.zip(weighted_densities) { let nwd = wd.shape()[0]; let ngrid = wd.len() / nwd; helmholtz_energy_density.push( - c.calculate_helmholtz_energy_density( + c.helmholtz_energy_density( temperature_dual, wd.into_shape((nwd, ngrid)).unwrap().view(), )? @@ -177,7 +174,7 @@ where temperature: Dual64, density: &Array, ) -> Array { - let lambda = self.dft.ideal_gas_model().ln_lambda3(temperature); + let lambda = self.dft.ln_lambda3(temperature); let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); for (i, rhoi) in density.outer_iter().enumerate() { phi += &rhoi.mapv(|rhoi| (lambda[i] + rhoi.ln() - 1.0) * rhoi); @@ -197,7 +194,6 @@ where let temperature_dual = Dual64::from(temperature).derivative(); let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions - .iter() .map(|c| c.weight_functions(temperature_dual)) .collect(); let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None); @@ -241,7 +237,6 @@ where let temperature_dual = Dual64::from(temperature).derivative(); let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions - .iter() .map(|c| c.weight_functions(temperature_dual)) .collect(); let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None); @@ -355,7 +350,6 @@ where // calculate temperature derivative of functional derivative let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions - .iter() .map(|c| c.weight_functions(Dual64::from(t).derivative())) .collect(); let convolver: Arc> = diff --git a/feos-dft/src/solver.rs b/feos-dft/src/solver.rs index 9f96c2449..ed0bd07c3 100644 --- a/feos-dft/src/solver.rs +++ b/feos-dft/src/solver.rs @@ -1,4 +1,7 @@ -use crate::{DFTProfile, HelmholtzEnergyFunctional, WeightFunction, WeightFunctionShape}; +use crate::{ + DFTProfile, FunctionalContribution, HelmholtzEnergyFunctional, WeightFunction, + WeightFunctionShape, +}; use feos_core::si::{Time, SECOND}; use feos_core::{log_iter, log_result, EosError, EosResult, Verbosity}; use nalgebra::{DMatrix, DVector}; @@ -559,8 +562,8 @@ where let temperature = self.temperature.to_reduced(); let contributions = self.dft.contributions(); let weighted_densities = self.convolver.weighted_densities(density); - let mut second_partial_derivatives = Vec::with_capacity(contributions.len()); - for (c, wd) in contributions.iter().zip(&weighted_densities) { + let mut second_partial_derivatives = Vec::new(); + for (c, wd) in contributions.zip(&weighted_densities) { let nwd = wd.shape()[0]; let ngrid = wd.len() / nwd; let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); diff --git a/src/association/dft.rs b/src/association/dft.rs index 5a8838a2c..9b58f984d 100644 --- a/src/association/dft.rs +++ b/src/association/dft.rs @@ -1,21 +1,15 @@ use super::*; use crate::hard_sphere::HardSphereProperties; use feos_core::EosResult; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, 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 + Copy + ScalarOperand, - P: HardSphereProperties, -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for Association

{ + fn weight_functions + Copy>(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; let r = p.hs_diameter(temperature) * 0.5; let [_, _, _, c3] = p.geometry_coefficients(temperature); @@ -42,7 +36,7 @@ where ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, @@ -105,15 +99,12 @@ where // auxiliary variables let n3i = n3.mapv(|n3| (-n3 + 1.0).recip()); - self.calculate_helmholtz_energy_density(temperature, &rho0, &n2, &n3i, &xi) + self._helmholtz_energy_density(temperature, &rho0, &n2, &n3i, &xi) } } impl Association

{ - pub fn calculate_helmholtz_energy_density< - N: DualNum + Copy + ScalarOperand, - S: Data, - >( + pub fn _helmholtz_energy_density + Copy + ScalarOperand, S: Data>( &self, temperature: N, rho0: &Array2, diff --git a/src/association/mod.rs b/src/association/mod.rs index d6f67492b..2426c967d 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -270,7 +270,7 @@ impl AssociationParameters { /// contribution and functional. pub struct Association

{ parameters: Arc

, - association_parameters: AssociationParameters, + association_parameters: Arc, max_iter: usize, tol: f64, force_cross_association: bool, @@ -285,7 +285,7 @@ impl Association

{ ) -> Self { Self { parameters: parameters.clone(), - association_parameters: association_parameters.clone(), + association_parameters: Arc::new(association_parameters.clone()), max_iter, tol, force_cross_association: false, @@ -352,7 +352,7 @@ impl Association

{ // association strength let [delta_ab, delta_cc] = - self.association_strength(state.temperature, &diameter, n2, n3i, D::one()); + self.association_strength(state.temperature, diameter, n2, n3i, D::one()); match ( a.sites_a.len() * a.sites_b.len(), diff --git a/src/eos.rs b/src/eos.rs index 150f6b7d9..9b23423d1 100644 --- a/src/eos.rs +++ b/src/eos.rs @@ -14,7 +14,7 @@ use feos_core::python::user_defined::PyResidual; use feos_core::si::*; use feos_core::*; use feos_derive::{Components, Residual}; -use ndarray::Array1; +use ndarray::{Array1, ScalarOperand}; use num_dual::DualNum; /// Collection of different [Residual] implementations. diff --git a/src/functional.rs b/src/functional.rs index a9b3c3c12..06cb2c6a8 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -1,19 +1,21 @@ #[cfg(feature = "gc_pcsaft")] -use crate::gc_pcsaft::GcPcSaftFunctional; -use crate::hard_sphere::FMTFunctional; +use crate::gc_pcsaft::{GcPcSaftFunctional, GcPcSaftFunctionalContribution}; +use crate::hard_sphere::{FMTContribution, FMTFunctional, HardSphereParameters}; #[cfg(feature = "pcsaft")] -use crate::pcsaft::PcSaftFunctional; +use crate::pcsaft::{PcSaftFunctional, PcSaftFunctionalContribution}; #[cfg(feature = "pets")] -use crate::pets::PetsFunctional; +use crate::pets::{PetsFunctional, PetsFunctionalContribution}; #[cfg(feature = "saftvrqmie")] -use crate::saftvrqmie::SaftVRQMieFunctional; +use crate::saftvrqmie::{SaftVRQMieFunctional, SaftVRQMieFunctionalContribution}; use feos_core::si::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}; +use ndarray::{Array1, Array2, ArrayView2, ScalarOperand}; +use num_dual::DualNum; use petgraph::graph::UnGraph; use petgraph::Graph; @@ -38,3 +40,16 @@ pub enum FunctionalVariant { #[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/dispersion.rs b/src/gc_pcsaft/dft/dispersion.rs index dfbeac2f4..e1a86914f 100644 --- a/src/gc_pcsaft/dft/dispersion.rs +++ b/src/gc_pcsaft/dft/dispersion.rs @@ -2,9 +2,7 @@ use super::parameter::GcPcSaftFunctionalParameters; use crate::gc_pcsaft::eos::dispersion::{A0, A1, A2, B0, B1, B2}; use crate::hard_sphere::HardSphereProperties; use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_6, PI}; @@ -24,10 +22,11 @@ impl AttractiveFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual - for AttractiveFunctional -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for AttractiveFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); @@ -37,7 +36,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, density: ArrayView2, diff --git a/src/gc_pcsaft/dft/hard_chain.rs b/src/gc_pcsaft/dft/hard_chain.rs index c3552821f..0eaeda1c6 100644 --- a/src/gc_pcsaft/dft/hard_chain.rs +++ b/src/gc_pcsaft/dft/hard_chain.rs @@ -1,9 +1,7 @@ use super::GcPcSaftFunctionalParameters; use crate::hard_sphere::HardSphereProperties; use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use petgraph::visit::EdgeRef; @@ -23,8 +21,11 @@ impl ChainFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual for ChainFunctional { - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for ChainFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); WeightFunctionInfo::new(p.component_index.clone(), true) @@ -46,7 +47,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for C ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 1fd181001..b3b600485 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -3,10 +3,13 @@ use crate::association::Association; use crate::hard_sphere::{FMTContribution, FMTVersion, HardSphereProperties, MonomerShape}; use feos_core::parameter::ParameterHetero; use feos_core::si::{MolarWeight, GRAM, MOL}; -use feos_core::Components; +use feos_core::{Components, EosResult}; +use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; -use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use ndarray::Array1; +use feos_dft::{ + FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, +}; +use ndarray::{Array1, ArrayView2, ScalarOperand}; use num_dual::DualNum; use petgraph::graph::UnGraph; use std::f64::consts::FRAC_PI_6; @@ -24,7 +27,6 @@ pub struct GcPcSaftFunctional { pub parameters: Arc, fmt_version: FMTVersion, options: GcPcSaftOptions, - contributions: Vec>, } impl GcPcSaftFunctional { @@ -41,36 +43,10 @@ impl GcPcSaftFunctional { fmt_version: FMTVersion, saft_options: GcPcSaftOptions, ) -> DFT { - let mut contributions: Vec> = Vec::with_capacity(4); - - // Hard sphere contribution - let hs = FMTContribution::new(¶meters, fmt_version); - contributions.push(Box::new(hs)); - - // Hard chains - let chain = ChainFunctional::new(¶meters); - contributions.push(Box::new(chain)); - - // Dispersion - let att = AttractiveFunctional::new(¶meters); - contributions.push(Box::new(att)); - - // Association - if !parameters.association.is_empty() { - let assoc = Association::new( - ¶meters, - ¶meters.association, - saft_options.max_iter_cross_assoc, - saft_options.tol_cross_assoc, - ); - contributions.push(Box::new(assoc)); - } - DFT(Self { parameters, fmt_version, options: saft_options, - contributions, }) } } @@ -91,6 +67,8 @@ impl Components for GcPcSaftFunctional { } impl HelmholtzEnergyFunctional for GcPcSaftFunctional { + type Contribution = GcPcSaftFunctionalContribution; + fn molecule_shape(&self) -> MoleculeShape { MoleculeShape::Heterosegmented(&self.parameters.component_index) } @@ -102,8 +80,33 @@ impl HelmholtzEnergyFunctional for GcPcSaftFunctional { / (FRAC_PI_6 * &p.m * p.sigma.mapv(|v| v.powi(3)) * moles_segments).sum() } - fn contributions(&self) -> &[Box] { - &self.contributions + fn contributions(&self) -> Box> { + let mut contributions = Vec::with_capacity(4); + + // Hard sphere contribution + let hs = FMTContribution::new(&self.parameters, self.fmt_version); + contributions.push(hs.into()); + + // Hard chains + let chain = ChainFunctional::new(&self.parameters); + contributions.push(chain.into()); + + // Dispersion + let att = AttractiveFunctional::new(&self.parameters); + contributions.push(att.into()); + + // Association + if !self.parameters.association.is_empty() { + let assoc = Association::new( + &self.parameters, + &self.parameters.association, + self.options.max_iter_cross_assoc, + self.options.tol_cross_assoc, + ); + contributions.push(Box::new(assoc).into()); + } + + Box::new(contributions.into_iter()) } fn molar_weight(&self) -> MolarWeight> { @@ -149,3 +152,11 @@ impl FluidParameters for GcPcSaftFunctional { &self.parameters.sigma } } + +#[derive(FunctionalContribution)] +pub enum GcPcSaftFunctionalContribution { + Fmt(FMTContribution), + Chain(ChainFunctional), + Attractive(AttractiveFunctional), + Association(Box>), +} diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index b3d7f93ae..32a7bb3f6 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -118,23 +118,23 @@ impl Residual for GcPcSaft { v.push(( "Hard Sphere".to_string(), - self.hard_sphere.helmholtz_energy(&state), + self.hard_sphere.helmholtz_energy(state), )); v.push(( "Hard Sphere".to_string(), - self.hard_chain.helmholtz_energy(&state), + self.hard_chain.helmholtz_energy(state), )); v.push(( "Dispersion".to_string(), - self.dispersion.helmholtz_energy(&state), + self.dispersion.helmholtz_energy(state), )); if let Some(dipole) = self.dipole.as_ref() { - v.push(("Dipole".to_string(), dipole.helmholtz_energy(&state))) + v.push(("Dipole".to_string(), dipole.helmholtz_energy(state))) } if let Some(association) = self.association.as_ref() { v.push(( "Association".to_string(), - association.helmholtz_energy(&state, &d), + association.helmholtz_energy(state, &d), )) } v diff --git a/src/gc_pcsaft/mod.rs b/src/gc_pcsaft/mod.rs index 7ddcd87a3..09345e552 100644 --- a/src/gc_pcsaft/mod.rs +++ b/src/gc_pcsaft/mod.rs @@ -12,7 +12,7 @@ pub(crate) mod eos; pub mod micelles; mod record; #[cfg(feature = "dft")] -pub use dft::{GcPcSaftFunctional, GcPcSaftFunctionalParameters}; +pub use dft::{GcPcSaftFunctional, GcPcSaftFunctionalContribution, GcPcSaftFunctionalParameters}; pub use eos::{GcPcSaft, GcPcSaftChemicalRecord, GcPcSaftEosParameters, GcPcSaftOptions}; pub use record::GcPcSaftRecord; diff --git a/src/hard_sphere/dft.rs b/src/hard_sphere/dft.rs index 886975815..2f4a2b7de 100644 --- a/src/hard_sphere/dft.rs +++ b/src/hard_sphere/dft.rs @@ -3,8 +3,8 @@ use feos_core::{Components, EosResult}; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; use feos_dft::{ - FunctionalContribution, FunctionalContributionDual, HelmholtzEnergyFunctional, MoleculeShape, - WeightFunction, WeightFunctionInfo, WeightFunctionShape, DFT, + FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunction, + WeightFunctionInfo, WeightFunctionShape, DFT, }; use ndarray::*; use num_dual::DualNum; @@ -85,10 +85,8 @@ impl

FMTContribution

{ } } -impl + Copy> FunctionalContributionDual - for FMTContribution

-{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for FMTContribution

{ + fn weight_functions + Copy>(&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()) { @@ -191,7 +189,7 @@ impl + Copy> FunctionalContributionDual } } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy>( &self, temperature: N, weighted_densities: ArrayView2, @@ -295,7 +293,7 @@ impl fmt::Display for FMTContribution

{ } } -struct HardSphereParameters { +pub struct HardSphereParameters { sigma: Array1, } @@ -312,7 +310,6 @@ impl HardSphereProperties for HardSphereParameters { /// [HelmholtzEnergyFunctional] for hard sphere systems. pub struct FMTFunctional { properties: Arc, - contributions: Vec>, version: FMTVersion, } @@ -321,11 +318,8 @@ impl FMTFunctional { let properties = Arc::new(HardSphereParameters { sigma: sigma.clone(), }); - let contributions: Vec> = - vec![Box::new(FMTContribution::new(&properties, version))]; DFT(Self { properties, - contributions, version, }) } @@ -346,8 +340,13 @@ impl Components for FMTFunctional { } impl HelmholtzEnergyFunctional for FMTFunctional { - fn contributions(&self) -> &[Box] { - &self.contributions + type Contribution = FMTContribution; + + fn contributions(&self) -> Box>> { + Box::new(std::iter::once(FMTContribution::new( + &self.properties, + self.version, + ))) } fn compute_max_density(&self, moles: &Array1) -> f64 { diff --git a/src/hard_sphere/mod.rs b/src/hard_sphere/mod.rs index 9aac82d99..eb2f0bd9f 100644 --- a/src/hard_sphere/mod.rs +++ b/src/hard_sphere/mod.rs @@ -10,7 +10,7 @@ use std::{borrow::Cow, sync::Arc}; #[cfg(feature = "dft")] mod dft; #[cfg(feature = "dft")] -pub use dft::{FMTContribution, FMTFunctional, FMTVersion}; +pub use dft::{FMTContribution, FMTFunctional, FMTVersion, HardSphereParameters}; /// Different monomer shapes for FMT and BMCSL. pub enum MonomerShape<'a, D> { diff --git a/src/ideal_gas/dippr.rs b/src/ideal_gas/dippr.rs index 27b70f3ac..e6356bd63 100644 --- a/src/ideal_gas/dippr.rs +++ b/src/ideal_gas/dippr.rs @@ -196,11 +196,9 @@ impl IdealGas for Dippr { }) .collect() } -} -impl fmt::Display for Dippr { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Ideal gas (DIPPR)") + fn ideal_gas_model(&self) -> String { + "Ideal gas (DIPPR)".into() } } diff --git a/src/ideal_gas/joback.rs b/src/ideal_gas/joback.rs index 523fc98bf..2d4498b93 100644 --- a/src/ideal_gas/joback.rs +++ b/src/ideal_gas/joback.rs @@ -144,11 +144,9 @@ impl IdealGas for Joback { (h - t * s) / (t * RGAS) + f }) } -} -impl fmt::Display for Joback { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Ideal gas (Joback)") + fn ideal_gas_model(&self) -> String { + "Ideal gas (Joback)".into() } } diff --git a/src/ideal_gas/mod.rs b/src/ideal_gas/mod.rs index a8a2ef9f6..3db7d22fc 100644 --- a/src/ideal_gas/mod.rs +++ b/src/ideal_gas/mod.rs @@ -1,13 +1,11 @@ //! Collection of ideal gas models. -use std::fmt; - #[cfg(feature = "python")] use feos_core::python::user_defined::PyIdealGas; use feos_core::{Components, IdealGas}; use feos_derive::{Components, IdealGas}; -use std::sync::Arc; use ndarray::Array1; use num_dual::DualNum; +use std::sync::Arc; mod dippr; mod joback; diff --git a/src/pcsaft/dft/dispersion.rs b/src/pcsaft/dft/dispersion.rs index 159199875..67fd85495 100644 --- a/src/pcsaft/dft/dispersion.rs +++ b/src/pcsaft/dft/dispersion.rs @@ -1,11 +1,9 @@ -use super::polar::calculate_helmholtz_energy_density_polar; +use super::polar::helmholtz_energy_density_polar; use super::PcSaftParameters; use crate::hard_sphere::HardSphereProperties; use crate::pcsaft::eos::dispersion::{A0, A1, A2, B0, B1, B2}; use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; @@ -40,18 +38,22 @@ fn att_weight_functions + Copy + ScalarOperand>( ) } -impl + Copy + ScalarOperand> FunctionalContributionDual - for AttractiveFunctional -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for AttractiveFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_DFT, temperature) } - fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_PDGT, temperature) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, density: ArrayView2, @@ -124,7 +126,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual }); // Helmholtz energy density - let phi_polar = calculate_helmholtz_energy_density_polar(p, temperature, density)?; + let phi_polar = helmholtz_energy_density_polar(p, temperature, density)?; Ok((-rho1mix * i1 * 2.0 - rho2mix * m_bar * c1 * i2) * PI + phi_polar) } } diff --git a/src/pcsaft/dft/hard_chain.rs b/src/pcsaft/dft/hard_chain.rs index 7119e2e40..3a8e2a632 100644 --- a/src/pcsaft/dft/hard_chain.rs +++ b/src/pcsaft/dft/hard_chain.rs @@ -1,9 +1,7 @@ use super::PcSaftParameters; use crate::hard_sphere::HardSphereProperties; use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use std::fmt; @@ -20,8 +18,11 @@ impl ChainFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual for ChainFunctional { - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for ChainFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let p = &self.parameters; let d = p.hs_diameter(temperature); WeightFunctionInfo::new(p.component_index().into_owned(), true) @@ -47,7 +48,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for C ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index 6729300de..5fa35414b 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -4,11 +4,15 @@ use crate::hard_sphere::{FMTContribution, FMTVersion}; use crate::pcsaft::eos::PcSaftOptions; use feos_core::parameter::Parameter; use feos_core::si::{MolarWeight, GRAM, MOL}; -use feos_core::Components; +use feos_core::{Components, EosResult}; +use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; -use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use ndarray::{Array1, Array2}; +use feos_dft::{ + FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, +}; +use ndarray::{Array1, Array2, ArrayView2, ScalarOperand}; +use num_dual::DualNum; use num_traits::One; use std::f64::consts::FRAC_PI_6; use std::sync::Arc; @@ -26,7 +30,6 @@ pub struct PcSaftFunctional { pub parameters: Arc, fmt_version: FMTVersion, options: PcSaftOptions, - contributions: Vec>, } impl PcSaftFunctional { @@ -43,53 +46,10 @@ impl PcSaftFunctional { fmt_version: FMTVersion, saft_options: PcSaftOptions, ) -> DFT { - let mut contributions: Vec> = Vec::with_capacity(4); - - if matches!( - fmt_version, - FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear - ) && parameters.m.len() == 1 - { - let fmt_assoc = PureFMTAssocFunctional::new(parameters.clone(), fmt_version); - contributions.push(Box::new(fmt_assoc)); - if parameters.m.iter().any(|&mi| mi > 1.0) { - let chain = PureChainFunctional::new(parameters.clone()); - contributions.push(Box::new(chain)); - } - let att = PureAttFunctional::new(parameters.clone()); - contributions.push(Box::new(att)); - } else { - // Hard sphere contribution - let hs = FMTContribution::new(¶meters, fmt_version); - contributions.push(Box::new(hs)); - - // Hard chains - if parameters.m.iter().any(|&mi| !mi.is_one()) { - let chain = ChainFunctional::new(parameters.clone()); - contributions.push(Box::new(chain)); - } - - // Dispersion - let att = AttractiveFunctional::new(parameters.clone()); - contributions.push(Box::new(att)); - - // Association - if !parameters.association.is_empty() { - let assoc = Association::new( - ¶meters, - ¶meters.association, - saft_options.max_iter_cross_assoc, - saft_options.tol_cross_assoc, - ); - contributions.push(Box::new(assoc)); - } - } - DFT(Self { parameters, fmt_version, options: saft_options, - contributions, }) } } @@ -110,14 +70,57 @@ impl Components for PcSaftFunctional { } impl HelmholtzEnergyFunctional for PcSaftFunctional { + type Contribution = PcSaftFunctionalContribution; + 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 contributions(&self) -> &[Box] { - &self.contributions + fn contributions(&self) -> Box> { + let mut contributions = Vec::with_capacity(4); + + if matches!( + self.fmt_version, + FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear + ) && self.parameters.m.len() == 1 + { + let fmt_assoc = PureFMTAssocFunctional::new(self.parameters.clone(), self.fmt_version); + contributions.push(fmt_assoc.into()); + if self.parameters.m.iter().any(|&mi| mi > 1.0) { + let chain = PureChainFunctional::new(self.parameters.clone()); + contributions.push(chain.into()); + } + let att = PureAttFunctional::new(self.parameters.clone()); + contributions.push(att.into()); + } else { + // Hard sphere contribution + let hs = FMTContribution::new(&self.parameters, self.fmt_version); + contributions.push(hs.into()); + + // Hard chains + if self.parameters.m.iter().any(|&mi| !mi.is_one()) { + let chain = ChainFunctional::new(self.parameters.clone()); + contributions.push(chain.into()); + } + + // Dispersion + let att = AttractiveFunctional::new(self.parameters.clone()); + contributions.push(att.into()); + + // Association + if !self.parameters.association.is_empty() { + let assoc = Association::new( + &self.parameters, + &self.parameters.association, + self.options.max_iter_cross_assoc, + self.options.tol_cross_assoc, + ); + contributions.push(assoc.into()); + } + } + Box::new(contributions.into_iter()) } fn molecule_shape(&self) -> MoleculeShape { @@ -149,3 +152,14 @@ impl PairPotential for PcSaftFunctional { }) } } + +#[derive(FunctionalContribution)] +pub enum PcSaftFunctionalContribution { + PureFMTAssoc(PureFMTAssocFunctional), + PureChain(PureChainFunctional), + PureAtt(PureAttFunctional), + Fmt(FMTContribution), + Chain(ChainFunctional), + Attractive(AttractiveFunctional), + Association(Association), +} diff --git a/src/pcsaft/dft/polar.rs b/src/pcsaft/dft/polar.rs index 16b0b4c10..c873d388a 100644 --- a/src/pcsaft/dft/polar.rs +++ b/src/pcsaft/dft/polar.rs @@ -8,7 +8,7 @@ use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; -pub(super) fn calculate_helmholtz_energy_density_polar + Copy + ScalarOperand>( +pub(super) fn helmholtz_energy_density_polar + Copy + ScalarOperand>( parameters: &PcSaftParameters, temperature: N, density: ArrayView2, diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index d8b0291a1..583e3a397 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -5,9 +5,7 @@ use crate::hard_sphere::{FMTVersion, HardSphereProperties}; 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_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; @@ -35,10 +33,11 @@ impl PureFMTAssocFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual - for PureFMTAssocFunctional -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for PureFMTAssocFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let r = self.parameters.hs_diameter(temperature) * 0.5; WeightFunctionInfo::new(arr1(&[0]), false).extend( vec![ @@ -57,7 +56,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, @@ -126,7 +125,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual }); let rho0 = (&n0 / p.m[0] * &xi).insert_axis(Axis(0)); - phi += &(self.association.calculate_helmholtz_energy_density( + phi += &(self.association._helmholtz_energy_density( temperature, &rho0, &n2, @@ -156,8 +155,11 @@ impl PureChainFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual for PureChainFunctional { - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for PureChainFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); WeightFunctionInfo::new(arr1(&[0]), true) .add( @@ -174,7 +176,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for P ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, _: N, weighted_densities: ArrayView2, @@ -208,8 +210,11 @@ impl PureAttFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual for PureAttFunctional { - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for PureAttFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); const PSI: f64 = 1.3862; // Homosegmented DFT (Sauer2017) WeightFunctionInfo::new(arr1(&[0]), false).add( @@ -218,7 +223,10 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for P ) } - fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); const PSI: f64 = 1.3286; // pDGT (Rehner2018) WeightFunctionInfo::new(arr1(&[0]), false).add( @@ -227,7 +235,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for P ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, diff --git a/src/pcsaft/eos/dispersion.rs b/src/pcsaft/eos/dispersion.rs index 5a82cb3ab..c0b052204 100644 --- a/src/pcsaft/eos/dispersion.rs +++ b/src/pcsaft/eos/dispersion.rs @@ -1,5 +1,4 @@ use super::PcSaftParameters; -use crate::hard_sphere::HardSphereProperties; use feos_core::StateHD; use ndarray::Array1; use num_dual::DualNum; @@ -67,7 +66,11 @@ pub struct Dispersion { } impl Dispersion { - pub fn helmholtz_energy + Copy>(&self, state: &StateHD, diameter: &Array1) -> D { + pub fn helmholtz_energy + Copy>( + &self, + state: &StateHD, + diameter: &Array1, + ) -> D { // auxiliary variables let n = self.parameters.m.len(); let p = &self.parameters; @@ -130,6 +133,7 @@ impl fmt::Display for Dispersion { #[cfg(test)] mod tests { use super::*; + use crate::hard_sphere::HardSphereProperties; use crate::pcsaft::parameters::utils::{ butane_parameters, propane_butane_parameters, propane_parameters, }; diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index 0854596ce..8a1211b33 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -4,7 +4,7 @@ use crate::hard_sphere::{HardSphere, HardSphereProperties}; use feos_core::parameter::Parameter; use feos_core::{si::*, StateHD}; use feos_core::{Components, EntropyScaling, EosError, EosResult, Residual, State}; -use ndarray::{Array1, ScalarOperand}; +use ndarray::Array1; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_6, PI}; use std::fmt; @@ -146,34 +146,34 @@ impl Residual for PcSaft { v.push(( "Hard Sphere".to_string(), - self.hard_sphere.helmholtz_energy(&state), + self.hard_sphere.helmholtz_energy(state), )); if let Some(hc) = self.hard_chain.as_ref() { - v.push(("Hard Sphere".to_string(), hc.helmholtz_energy(&state))) + v.push(("Hard Chain".to_string(), hc.helmholtz_energy(state))) } v.push(( "Dispersion".to_string(), - self.dispersion.helmholtz_energy(&state, &d), + self.dispersion.helmholtz_energy(state, &d), )); if let Some(dipole) = self.dipole.as_ref() { - v.push(("Dipole".to_string(), dipole.helmholtz_energy(&state, &d))) + v.push(("Dipole".to_string(), dipole.helmholtz_energy(state, &d))) } if let Some(quadrupole) = self.quadrupole.as_ref() { v.push(( "Quadrupole".to_string(), - quadrupole.helmholtz_energy(&state, &d), + quadrupole.helmholtz_energy(state, &d), )) } if let Some(dipole_quadrupole) = self.dipole_quadrupole.as_ref() { v.push(( "Dipole-Quadrupole".to_string(), - dipole_quadrupole.helmholtz_energy(&state, &d), + dipole_quadrupole.helmholtz_energy(state, &d), )) } if let Some(association) = self.association.as_ref() { v.push(( "Association".to_string(), - association.helmholtz_energy(&state, &d), + association.helmholtz_energy(state, &d), )) } v diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs index 7e32eaa10..105729cf5 100644 --- a/src/pcsaft/eos/polar.rs +++ b/src/pcsaft/eos/polar.rs @@ -1,5 +1,4 @@ use super::PcSaftParameters; -use crate::hard_sphere::HardSphereProperties; use feos_core::StateHD; use ndarray::prelude::*; use num_dual::DualNum; @@ -167,7 +166,11 @@ pub struct Dipole { impl Dipole { #[inline] - pub fn helmholtz_energy + Copy>(&self, state: &StateHD, diameter: &Array1) -> D { + pub fn helmholtz_energy + Copy>( + &self, + state: &StateHD, + diameter: &Array1, + ) -> D { let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Dipole); let p = &self.parameters; @@ -246,7 +249,11 @@ pub struct Quadrupole { impl Quadrupole { #[inline] - pub fn helmholtz_energy + Copy>(&self, state: &StateHD, diameter: &Array1) -> D { + pub fn helmholtz_energy + Copy>( + &self, + state: &StateHD, + diameter: &Array1, + ) -> D { let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Quadrupole); let p = &self.parameters; @@ -335,7 +342,11 @@ pub struct DipoleQuadrupole { impl DipoleQuadrupole { #[inline] - pub fn helmholtz_energy + Copy>(&self, state: &StateHD, diameter: &Array1) -> D { + pub fn helmholtz_energy + Copy>( + &self, + state: &StateHD, + diameter: &Array1, + ) -> D { let p = &self.parameters; let t_inv = state.temperature.inv(); @@ -444,6 +455,7 @@ impl fmt::Display for DipoleQuadrupole { #[cfg(test)] mod tests { use super::*; + use crate::hard_sphere::HardSphereProperties; use crate::pcsaft::eos::dispersion::Dispersion; use crate::pcsaft::parameters::utils::{ carbon_dioxide_parameters, dme_co2_parameters, dme_parameters, diff --git a/src/pcsaft/mod.rs b/src/pcsaft/mod.rs index 2c36ea4ac..5d3744b4f 100644 --- a/src/pcsaft/mod.rs +++ b/src/pcsaft/mod.rs @@ -10,7 +10,7 @@ mod eos; pub(crate) mod parameters; #[cfg(feature = "dft")] -pub use dft::PcSaftFunctional; +pub use dft::{PcSaftFunctional, PcSaftFunctionalContribution}; pub use eos::{DQVariants, PcSaft, PcSaftOptions}; pub use parameters::{PcSaftBinaryRecord, PcSaftParameters, PcSaftRecord}; diff --git a/src/pets/dft/dispersion.rs b/src/pets/dft/dispersion.rs index e1e2072ad..27d3fbd50 100644 --- a/src/pets/dft/dispersion.rs +++ b/src/pets/dft/dispersion.rs @@ -2,9 +2,7 @@ use crate::hard_sphere::HardSphereProperties; use crate::pets::eos::dispersion::{A, B}; use crate::pets::parameters::PetsParameters; use feos_core::EosError; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use std::f64::consts::{FRAC_PI_3, PI}; @@ -39,18 +37,22 @@ fn att_weight_functions + Copy + ScalarOperand>( ) } -impl + Copy + ScalarOperand> FunctionalContributionDual - for AttractiveFunctional -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for AttractiveFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_DFT, temperature) } - fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_PDGT, temperature) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, density: ArrayView2, diff --git a/src/pets/dft/mod.rs b/src/pets/dft/mod.rs index 94002dff6..37040766e 100644 --- a/src/pets/dft/mod.rs +++ b/src/pets/dft/mod.rs @@ -4,11 +4,14 @@ use crate::hard_sphere::{FMTContribution, FMTVersion}; use dispersion::AttractiveFunctional; use feos_core::parameter::Parameter; use feos_core::si::{MolarWeight, GRAM, MOL}; -use feos_core::Components; +use feos_core::{Components, EosResult}; +use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; -use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use ndarray::{Array1, Array2}; +use feos_dft::{ + FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, +}; +use ndarray::{Array1, Array2, ArrayView2, ScalarOperand}; use num_dual::DualNum; use pure_pets_functional::*; use std::f64::consts::FRAC_PI_6; @@ -23,7 +26,6 @@ pub struct PetsFunctional { pub parameters: Arc, fmt_version: FMTVersion, options: PetsOptions, - contributions: Vec>, } impl PetsFunctional { @@ -47,36 +49,10 @@ impl PetsFunctional { fmt_version: FMTVersion, pets_options: PetsOptions, ) -> DFT { - let mut contributions: Vec> = Vec::with_capacity(2); - - if matches!( - fmt_version, - FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear - ) && parameters.sigma.len() == 1 - // Pure substance or mixture - { - // Hard-sphere contribution pure substance - let fmt = PureFMTFunctional::new(parameters.clone(), fmt_version); - contributions.push(Box::new(fmt)); - - // Dispersion contribution pure substance - let att = PureAttFunctional::new(parameters.clone()); - contributions.push(Box::new(att)); - } else { - // Hard-sphere contribution mixtures - let hs = FMTContribution::new(¶meters, fmt_version); - contributions.push(Box::new(hs)); - - // Dispersion contribution mixtures - let att = AttractiveFunctional::new(parameters.clone()); - contributions.push(Box::new(att)); - } - DFT(Self { parameters, fmt_version, options: pets_options, - contributions, }) } } @@ -97,6 +73,8 @@ impl Components for PetsFunctional { } impl HelmholtzEnergyFunctional for PetsFunctional { + type Contribution = PetsFunctionalContribution; + fn molecule_shape(&self) -> MoleculeShape { MoleculeShape::Spherical(self.parameters.sigma.len()) } @@ -106,8 +84,33 @@ impl HelmholtzEnergyFunctional for PetsFunctional { / (FRAC_PI_6 * self.parameters.sigma.mapv(|v| v.powi(3)) * moles).sum() } - fn contributions(&self) -> &[Box] { - &self.contributions + fn contributions(&self) -> Box<(dyn Iterator)> { + let mut contributions = Vec::with_capacity(2); + + if matches!( + self.fmt_version, + FMTVersion::WhiteBear | FMTVersion::AntiSymWhiteBear + ) && self.parameters.sigma.len() == 1 + // Pure substance or mixture + { + // Hard-sphere contribution pure substance + let fmt = PureFMTFunctional::new(self.parameters.clone(), self.fmt_version); + contributions.push(fmt.into()); + + // Dispersion contribution pure substance + let att = PureAttFunctional::new(self.parameters.clone()); + contributions.push(att.into()); + } else { + // Hard-sphere contribution mixtures + let hs = FMTContribution::new(&self.parameters, self.fmt_version); + contributions.push(hs.into()); + + // Dispersion contribution mixtures + let att = AttractiveFunctional::new(self.parameters.clone()); + contributions.push(att.into()); + } + + Box::new(contributions.into_iter()) } fn molar_weight(&self) -> MolarWeight> { @@ -140,3 +143,11 @@ impl PairPotential for PetsFunctional { }) } } + +#[derive(FunctionalContribution)] +pub enum PetsFunctionalContribution { + PureFMT(PureFMTFunctional), + PureAtt(PureAttFunctional), + Fmt(FMTContribution), + Attractive(AttractiveFunctional), +} diff --git a/src/pets/dft/pure_pets_functional.rs b/src/pets/dft/pure_pets_functional.rs index 85ec7c910..3e1e99d2d 100644 --- a/src/pets/dft/pure_pets_functional.rs +++ b/src/pets/dft/pure_pets_functional.rs @@ -2,9 +2,7 @@ use crate::hard_sphere::{FMTVersion, HardSphereProperties}; use crate::pets::eos::dispersion::{A, B}; use crate::pets::parameters::PetsParameters; use feos_core::{EosError, EosResult}; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::*; use std::f64::consts::{FRAC_PI_6, PI}; @@ -29,8 +27,11 @@ impl PureFMTFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual for PureFMTFunctional { - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for PureFMTFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let r = self.parameters.hs_diameter(temperature) * 0.5; WeightFunctionInfo::new(arr1(&[0]), false).extend( vec![ @@ -49,7 +50,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for P ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, @@ -125,8 +126,11 @@ impl PureAttFunctional { } } -impl + Copy + ScalarOperand> FunctionalContributionDual for PureAttFunctional { - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for PureAttFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); const PSI: f64 = 1.21; // Homosegmented DFT (Heier2018) WeightFunctionInfo::new(arr1(&[0]), false).add( @@ -135,7 +139,10 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for P ) } - fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { let d = self.parameters.hs_diameter(temperature); const PSI: f64 = 1.21; // pDGT (not yet determined) WeightFunctionInfo::new(arr1(&[0]), false).add( @@ -144,7 +151,7 @@ impl + Copy + ScalarOperand> FunctionalContributionDual for P ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, diff --git a/src/pets/mod.rs b/src/pets/mod.rs index c79c12580..5f496b9b4 100644 --- a/src/pets/mod.rs +++ b/src/pets/mod.rs @@ -16,7 +16,7 @@ mod eos; mod parameters; #[cfg(feature = "dft")] -pub use dft::PetsFunctional; +pub use dft::{PetsFunctional, PetsFunctionalContribution}; pub use eos::{Pets, PetsOptions}; pub use parameters::{PetsBinaryRecord, PetsParameters, PetsRecord}; diff --git a/src/saftvrqmie/dft/dispersion.rs b/src/saftvrqmie/dft/dispersion.rs index 4b21d0538..7acb07609 100644 --- a/src/saftvrqmie/dft/dispersion.rs +++ b/src/saftvrqmie/dft/dispersion.rs @@ -1,9 +1,7 @@ use crate::saftvrqmie::eos::dispersion::{dispersion_energy_density, Alpha}; use crate::saftvrqmie::parameters::SaftVRQMieParameters; use feos_core::EosResult; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use std::fmt; @@ -38,18 +36,22 @@ fn att_weight_functions + Copy + ScalarOperand>( ) } -impl + Copy + ScalarOperand> FunctionalContributionDual - for AttractiveFunctional -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for AttractiveFunctional { + fn weight_functions + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_DFT, temperature) } - fn weight_functions_pdgt(&self, temperature: N) -> WeightFunctionInfo { + fn weight_functions_pdgt + Copy + ScalarOperand>( + &self, + temperature: N, + ) -> WeightFunctionInfo { att_weight_functions(&self.parameters, PSI_PDGT, temperature) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, density: ArrayView2, diff --git a/src/saftvrqmie/dft/mod.rs b/src/saftvrqmie/dft/mod.rs index 3f57290b4..f972601b1 100644 --- a/src/saftvrqmie/dft/mod.rs +++ b/src/saftvrqmie/dft/mod.rs @@ -4,11 +4,14 @@ use crate::saftvrqmie::parameters::SaftVRQMieParameters; use dispersion::AttractiveFunctional; use feos_core::parameter::Parameter; use feos_core::si::{MolarWeight, GRAM, MOL}; -use feos_core::Components; +use feos_core::{Components, EosResult}; +use feos_derive::FunctionalContribution; use feos_dft::adsorption::FluidParameters; use feos_dft::solvation::PairPotential; -use feos_dft::{FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, DFT}; -use ndarray::{Array, Array1, Array2}; +use feos_dft::{ + FunctionalContribution, HelmholtzEnergyFunctional, MoleculeShape, WeightFunctionInfo, DFT, +}; +use ndarray::{Array, Array1, Array2, ArrayView2, ScalarOperand}; use non_additive_hs::NonAddHardSphereFunctional; use num_dual::DualNum; use std::f64::consts::FRAC_PI_6; @@ -22,7 +25,6 @@ pub struct SaftVRQMieFunctional { pub parameters: Arc, fmt_version: FMTVersion, options: SaftVRQMieOptions, - contributions: Vec>, } impl SaftVRQMieFunctional { @@ -43,27 +45,10 @@ impl SaftVRQMieFunctional { fmt_version: FMTVersion, saft_options: SaftVRQMieOptions, ) -> DFT { - let mut contributions: Vec> = Vec::with_capacity(3); - - // Hard sphere contribution - let hs = FMTContribution::new(¶meters, fmt_version); - contributions.push(Box::new(hs)); - - // Non-additive hard-sphere contribution - if saft_options.inc_nonadd_term { - let non_add_hs = NonAddHardSphereFunctional::new(parameters.clone()); - contributions.push(Box::new(non_add_hs)); - } - - // Dispersion - let att = AttractiveFunctional::new(parameters.clone()); - contributions.push(Box::new(att)); - DFT(Self { parameters, fmt_version, options: saft_options, - contributions, }) } } @@ -84,14 +69,32 @@ impl Components for SaftVRQMieFunctional { } impl HelmholtzEnergyFunctional for SaftVRQMieFunctional { + type Contribution = SaftVRQMieFunctionalContribution; + 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 contributions(&self) -> &[Box] { - &self.contributions + fn contributions(&self) -> Box<(dyn Iterator)> { + let mut contributions = Vec::with_capacity(3); + + // Hard sphere contribution + let hs = FMTContribution::new(&self.parameters, self.fmt_version); + contributions.push(hs.into()); + + // Non-additive hard-sphere contribution + if self.options.inc_nonadd_term { + let non_add_hs = NonAddHardSphereFunctional::new(self.parameters.clone()); + contributions.push(non_add_hs.into()); + } + + // Dispersion + let att = AttractiveFunctional::new(self.parameters.clone()); + contributions.push(att.into()); + + Box::new(contributions.into_iter()) } fn molar_weight(&self) -> MolarWeight> { @@ -130,3 +133,10 @@ impl PairPotential for SaftVRQMieFunctional { }) } } + +#[derive(FunctionalContribution)] +pub enum SaftVRQMieFunctionalContribution { + Fmt(FMTContribution), + NonAddHardSphere(NonAddHardSphereFunctional), + Attractive(AttractiveFunctional), +} diff --git a/src/saftvrqmie/dft/non_additive_hs.rs b/src/saftvrqmie/dft/non_additive_hs.rs index 04b5e7885..b9c3a81a7 100644 --- a/src/saftvrqmie/dft/non_additive_hs.rs +++ b/src/saftvrqmie/dft/non_additive_hs.rs @@ -1,8 +1,6 @@ use crate::saftvrqmie::parameters::SaftVRQMieParameters; use feos_core::EosResult; -use feos_dft::{ - FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape, -}; +use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use ndarray::*; use num_dual::DualNum; use std::f64::consts::PI; @@ -22,11 +20,11 @@ impl NonAddHardSphereFunctional { } } -impl FunctionalContributionDual for NonAddHardSphereFunctional -where - N: DualNum + Copy + ScalarOperand, -{ - fn weight_functions(&self, temperature: N) -> WeightFunctionInfo { +impl FunctionalContribution for NonAddHardSphereFunctional { + fn weight_functions + Copy + ScalarOperand>( + &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) @@ -52,7 +50,7 @@ where ) } - fn calculate_helmholtz_energy_density( + fn helmholtz_energy_density + Copy + ScalarOperand>( &self, temperature: N, weighted_densities: ArrayView2, diff --git a/src/saftvrqmie/eos/dispersion.rs b/src/saftvrqmie/eos/dispersion.rs index 6690f4286..35bc012c1 100644 --- a/src/saftvrqmie/eos/dispersion.rs +++ b/src/saftvrqmie/eos/dispersion.rs @@ -139,17 +139,17 @@ impl Dispersion { rho_s += rho[i] * p.m[i]; } // packing fractions - let zeta = zeta_saft_vrq_mie(&p.m, &x_s, &d_hs_ij, rho_s); - let zeta_bar = zeta_saft_vrq_mie(&p.m, &x_s, &s_eff_ij, rho_s); + let zeta = zeta_saft_vrq_mie(&p.m, &x_s, d_hs_ij, rho_s); + let zeta_bar = zeta_saft_vrq_mie(&p.m, &x_s, s_eff_ij, rho_s); // alphas .... - let alpha = Alpha::new(p, &s_eff_ij, &epsilon_k_eff_ij, state.temperature); + let alpha = Alpha::new(p, s_eff_ij, epsilon_k_eff_ij, state.temperature); - let a1 = first_order_perturbation(p, &x_s, zeta, rho_s, &d_hs_ij, &s_eff_ij, &dq_ij); + let a1 = first_order_perturbation(p, &x_s, zeta, rho_s, d_hs_ij, s_eff_ij, dq_ij); let a2 = second_order_perturbation( - p, &alpha, &x_s, zeta, zeta_bar, rho_s, &d_hs_ij, &s_eff_ij, &dq_ij, + p, &alpha, &x_s, zeta, zeta_bar, rho_s, d_hs_ij, s_eff_ij, dq_ij, ); - let a3 = third_order_perturbation(p, &alpha, &x_s, zeta_bar, &epsilon_k_eff_ij); + let a3 = third_order_perturbation(p, &alpha, &x_s, zeta_bar, epsilon_k_eff_ij); let mut n_s = D::zero(); for i in 0..n { diff --git a/src/saftvrqmie/eos/mod.rs b/src/saftvrqmie/eos/mod.rs index cc619764a..849517cfa 100644 --- a/src/saftvrqmie/eos/mod.rs +++ b/src/saftvrqmie/eos/mod.rs @@ -168,16 +168,16 @@ impl Residual for SaftVRQMie { v.push(( "Hard Sphere".to_string(), - self.hard_sphere.helmholtz_energy(&state, &properties), + self.hard_sphere.helmholtz_energy(state, &properties), )); v.push(( "Dispersion".to_string(), - self.dispersion.helmholtz_energy(&state, &properties), + self.dispersion.helmholtz_energy(state, &properties), )); if let Some(non_additive_hard_sphere) = self.non_additive_hard_sphere.as_ref() { v.push(( "Non additive Hard Sphere".to_string(), - non_additive_hard_sphere.helmholtz_energy(&state, &properties), + non_additive_hard_sphere.helmholtz_energy(state, &properties), )) } v diff --git a/src/saftvrqmie/eos/non_additive_hs.rs b/src/saftvrqmie/eos/non_additive_hs.rs index 2579d3d6d..806fbd605 100644 --- a/src/saftvrqmie/eos/non_additive_hs.rs +++ b/src/saftvrqmie/eos/non_additive_hs.rs @@ -28,7 +28,7 @@ impl NonAddHardSphere { Array2::from_shape_fn((n, n), |(i, j)| (d_hs_ij[[i, i]] + d_hs_ij[[j, j]]) * 0.5); let n_s = Array1::from_shape_fn(n, |i| state.moles[i] * p.m[i]).sum(); - n_s * reduced_non_additive_hs_energy(p, &d_hs_ij, &d_hs_add_ij, &state.partial_density) + n_s * reduced_non_additive_hs_energy(p, d_hs_ij, &d_hs_add_ij, &state.partial_density) } } diff --git a/src/saftvrqmie/mod.rs b/src/saftvrqmie/mod.rs index e56adc16c..158094933 100644 --- a/src/saftvrqmie/mod.rs +++ b/src/saftvrqmie/mod.rs @@ -15,7 +15,7 @@ mod eos; mod parameters; #[cfg(feature = "dft")] -pub use dft::SaftVRQMieFunctional; +pub use dft::{SaftVRQMieFunctional, SaftVRQMieFunctionalContribution}; pub use eos::{FeynmanHibbsOrder, SaftVRQMie, SaftVRQMieOptions}; pub use parameters::{SaftVRQMieBinaryRecord, SaftVRQMieParameters, SaftVRQMieRecord}; diff --git a/src/uvtheory/eos/bh/attractive_perturbation.rs b/src/uvtheory/eos/bh/attractive_perturbation.rs index fbd7b0731..513678e1b 100644 --- a/src/uvtheory/eos/bh/attractive_perturbation.rs +++ b/src/uvtheory/eos/bh/attractive_perturbation.rs @@ -87,7 +87,11 @@ fn delta_b12u>(t_x: D, mean_field_constant_x: D, weighted_sigma3 -mean_field_constant_x / t_x * 2.0 * PI * weighted_sigma3_ij } -fn residual_virial_coefficient + Copy>(p: &UVTheoryParameters, x: &Array1, t: D) -> D { +fn residual_virial_coefficient + Copy>( + p: &UVTheoryParameters, + x: &Array1, + t: D, +) -> D { let mut delta_b2bar = D::zero(); for i in 0..p.ncomponents { let xi = x[i]; diff --git a/src/uvtheory/eos/bh/reference_perturbation.rs b/src/uvtheory/eos/bh/reference_perturbation.rs index ec1b6bc0b..31f1d9c30 100644 --- a/src/uvtheory/eos/bh/reference_perturbation.rs +++ b/src/uvtheory/eos/bh/reference_perturbation.rs @@ -1,6 +1,4 @@ -use super::hard_sphere::{ - diameter_bh, packing_fraction, packing_fraction_a, packing_fraction_b, -}; +use super::hard_sphere::{diameter_bh, packing_fraction, packing_fraction_a, packing_fraction_b}; use crate::uvtheory::parameters::*; use feos_core::StateHD; use num_dual::DualNum; diff --git a/src/uvtheory/eos/wca/attractive_perturbation.rs b/src/uvtheory/eos/wca/attractive_perturbation.rs index 997710cf9..d03983e86 100644 --- a/src/uvtheory/eos/wca/attractive_perturbation.rs +++ b/src/uvtheory/eos/wca/attractive_perturbation.rs @@ -123,7 +123,11 @@ fn delta_b12u + Copy>( * weighted_sigma3_ij } -fn residual_virial_coefficient + Copy>(p: &UVTheoryParameters, x: &Array1, t: D) -> D { +fn residual_virial_coefficient + Copy>( + p: &UVTheoryParameters, + x: &Array1, + t: D, +) -> D { let mut delta_b2bar = D::zero(); for i in 0..p.ncomponents { diff --git a/src/uvtheory/eos/wca/attractive_perturbation_uvb3.rs b/src/uvtheory/eos/wca/attractive_perturbation_uvb3.rs index 8aaa2cc7f..6f5a98f34 100644 --- a/src/uvtheory/eos/wca/attractive_perturbation_uvb3.rs +++ b/src/uvtheory/eos/wca/attractive_perturbation_uvb3.rs @@ -144,7 +144,11 @@ fn delta_b12u + Copy>( * weighted_sigma3_ij } -fn residual_virial_coefficient + Copy>(p: &UVTheoryParameters, x: &Array1, t: D) -> D { +fn residual_virial_coefficient + Copy>( + p: &UVTheoryParameters, + x: &Array1, + t: D, +) -> D { let mut delta_b2bar = D::zero(); for i in 0..p.ncomponents { diff --git a/src/uvtheory/parameters.rs b/src/uvtheory/parameters.rs index 0dd11e5dd..d0509a9bf 100644 --- a/src/uvtheory/parameters.rs +++ b/src/uvtheory/parameters.rs @@ -198,7 +198,12 @@ impl Parameter for UVTheoryParameters { }) } - fn records(&self) -> (&[PureRecord], Option<&Array2>) { + fn records( + &self, + ) -> ( + &[PureRecord], + Option<&Array2>, + ) { (&self.pure_records, self.binary_records.as_ref()) } }