diff --git a/docs/api/dft.md b/docs/api/dft.md index bcd65b5d6..b93b56d1b 100644 --- a/docs/api/dft.md +++ b/docs/api/dft.md @@ -28,11 +28,13 @@ Implementations of Helmholtz energy functionals for DFT. :toctree: generated/ State + StateVec PhaseEquilibrium PhaseDiagram Contributions Verbosity FMTVersion + DFTSolver ``` ## Interfaces @@ -54,6 +56,7 @@ Implementations of Helmholtz energy functionals for DFT. ExternalPotential Geometry Pore1D + Pore2D Pore3D Adsorption1D Adsorption3D diff --git a/docs/api/gc_pcsaft.md b/docs/api/gc_pcsaft.md index 6a6d68e45..ce72457df 100644 --- a/docs/api/gc_pcsaft.md +++ b/docs/api/gc_pcsaft.md @@ -1,6 +1,6 @@ # `feos.gc_pcsaft` -Utilities to build `GcPcSaftParameters`. To learn more about ways to build parameters from files or within Python, see [this example](/examples/eos/pcsaft/pcsaft_working_with_parameters). +Utilities to build `GcPcSaftParameters`. To learn more about ways to build parameters from files or within Python, see [this example](/tutorials/eos/pcsaft/pcsaft_working_with_parameters). ## Example @@ -20,10 +20,10 @@ from feos.gc_pcsaft import GcPcSaftParameters Identifier IdentifierOption ChemicalRecord - AssociationRecord + SmartsRecord + GcPcSaftRecord SegmentRecord BinarySegmentRecord - GcPcSaftRecord GcPcSaftEosParameters GcPcSaftFunctionalParameters ``` \ No newline at end of file diff --git a/docs/api/pcsaft.md b/docs/api/pcsaft.md index 63d6784fa..70166dd81 100644 --- a/docs/api/pcsaft.md +++ b/docs/api/pcsaft.md @@ -1,6 +1,6 @@ # `feos.pcsaft` -Utilities to build `PcSaftParameters`. To learn more about ways to build parameters from files or within Python, see [this example](/examples/eos/pcsaft/pcsaft_working_with_parameters). +Utilities to build `PcSaftParameters`. To learn more about ways to build parameters from files or within Python, see [this example](/tutorials/eos/pcsaft/pcsaft_working_with_parameters). ## Example @@ -19,12 +19,15 @@ parameters = PcSaftParameters.from_json(['methane', 'ethane'], 'parameters.json' :toctree: generated/ Identifier + IdentifierOption ChemicalRecord + SmartsRecord + DQVariants + PcSaftRecord + PcSaftBinaryRecord PureRecord SegmentRecord BinaryRecord BinarySegmentRecord - DQVariants - PcSaftRecord PcSaftParameters ``` \ No newline at end of file diff --git a/docs/api/pets.md b/docs/api/pets.md index 87c7e0e40..c3272f532 100644 --- a/docs/api/pets.md +++ b/docs/api/pets.md @@ -12,6 +12,7 @@ :toctree: generated/ Identifier + IdentifierOption ChemicalRecord PureRecord BinaryRecord diff --git a/docs/api/saftvrmie.md b/docs/api/saftvrmie.md index e9b5be7e0..77bb1d5c8 100644 --- a/docs/api/saftvrmie.md +++ b/docs/api/saftvrmie.md @@ -20,9 +20,10 @@ parameters = SaftVRMieParameters.from_json(['ethane', 'methane'], path) :toctree: generated/ Identifier - ChemicalRecord + IdentifierOption PureRecord BinaryRecord SaftVRMieRecord + SaftVRMieBinaryRecord SaftVRMieParameters ``` diff --git a/docs/api/saftvrqmie.md b/docs/api/saftvrqmie.md index e0ae90ca2..09454be22 100644 --- a/docs/api/saftvrqmie.md +++ b/docs/api/saftvrqmie.md @@ -20,6 +20,7 @@ parameters = SaftVRQMieParameters.from_json(['hydrogen', 'neon'], 'parameters.js FeynmanHibbsOrder Identifier + IdentifierOption PureRecord BinaryRecord SaftVRQMieRecord diff --git a/docs/api/si.md b/docs/api/si.md index fe1a5957f..01bce2b67 100644 --- a/docs/api/si.md +++ b/docs/api/si.md @@ -1,7 +1,7 @@ # `feos.si` This module contains data types for dimensioned quantities, i.e. floating point values multiplied with units. -If you want to learn more about this module, take a look at [this notebook](/examples/utility/working_with_units). +If you want to learn more about this module, take a look at [this notebook](/tutorials/utility/core_working_with_units). ## Example usage @@ -11,7 +11,7 @@ from feos.si import KELVIN, RGAS, METER, MOL, BAR p_ig = 25.0 * MOL / METER**3 * RGAS * 315.5 * KELVIN print('Ideal gas pressure: {:.2f} bar'.format(p_ig / BAR)) ``` -```terminal +``` Ideal gas pressure: 0.66 bar ``` diff --git a/docs/api/uvtheory.md b/docs/api/uvtheory.md index 9d4b0f05b..33775ef61 100644 --- a/docs/api/uvtheory.md +++ b/docs/api/uvtheory.md @@ -19,6 +19,7 @@ parameters = UVTheoryParameters.from_json(['methane', 'ethane'], 'parameters.jso :toctree: generated/ Identifier + IdentifierOption ChemicalRecord PureRecord BinaryRecord diff --git a/docs/recipes/index.md b/docs/recipes/index.md index 135edf2c5..28f64f143 100644 --- a/docs/recipes/index.md +++ b/docs/recipes/index.md @@ -1,7 +1,7 @@ # Recipes This section contains short code snippets for specific, commonly used tasks. -If you are looking for examples with explanations, see the [examples](/examples/index). +If you are looking for tutorials with explanations, see the [tutorials](/tutorials/index). ## Plots diff --git a/src/association/dft.rs b/src/association/dft.rs index 9b58f984d..0a5fb69be 100644 --- a/src/association/dft.rs +++ b/src/association/dft.rs @@ -1,5 +1,4 @@ use super::*; -use crate::hard_sphere::HardSphereProperties; use feos_core::EosResult; use feos_dft::{FunctionalContribution, WeightFunction, WeightFunctionInfo, WeightFunctionShape}; use num_dual::DualNum; @@ -8,7 +7,10 @@ use std::ops::MulAssign; pub const N0_CUTOFF: f64 = 1e-9; -impl FunctionalContribution for Association

{ +impl FunctionalContribution for Association

+where + P::Record: Sync + Send, +{ fn weight_functions + Copy>(&self, temperature: N) -> WeightFunctionInfo { let p = &self.parameters; let r = p.hs_diameter(temperature) * 0.5; @@ -103,7 +105,7 @@ impl FunctionalContribution for Associati } } -impl Association

{ +impl Association

{ pub fn _helmholtz_energy_density + Copy + ScalarOperand, S: Data>( &self, temperature: N, @@ -199,7 +201,9 @@ impl Association

{ let dj = diameter[j]; let k = n2 * n3i * (di * dj / (di + dj)); let delta = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) - * ((temperature.recip() * a.epsilon_k_ab[(0, 0)]).exp_m1() * a.sigma3_kappa_ab[(0, 0)]); + * self + .parameters + .association_strength(temperature, 0, 0, a.parameters_ab[(0, 0)]); // no cross association, two association sites let aux = &delta * (&rhob - &rhoa) + 1.0; @@ -233,7 +237,9 @@ impl Association

{ let di = diameter[i]; let k = n2 * n3i * (di * 0.5); let delta = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) - * ((temperature.recip() * a.epsilon_k_cc[(0, 0)]).exp_m1() * a.sigma3_kappa_cc[(0, 0)]); + * self + .parameters + .association_strength(temperature, 0, 0, a.parameters_cc[(0, 0)]); // no cross association, two association sites let xc = ((delta * 4.0 * &rhoc + 1.0).map(N::sqrt) + 1.0).map(N::recip) * 2.0; diff --git a/src/association/mod.rs b/src/association/mod.rs index dfae2a8d3..272d56dc5 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -13,43 +13,31 @@ use std::sync::Arc; #[cfg(feature = "dft")] mod dft; -#[cfg(feature = "python")] -mod python; -#[cfg(feature = "python")] -pub use python::PyAssociationRecord; #[derive(Clone, Copy, Debug)] -struct AssociationSite { +struct AssociationSite { assoc_comp: usize, site_index: usize, n: f64, - kappa_ab: f64, - epsilon_k_ab: f64, + parameters: A, } -impl AssociationSite { - fn new(assoc_comp: usize, site_index: usize, n: f64, kappa_ab: f64, epsilon_k_ab: f64) -> Self { +impl AssociationSite { + fn new(assoc_comp: usize, site_index: usize, n: f64, parameters: A) -> Self { Self { assoc_comp, site_index, n, - kappa_ab, - epsilon_k_ab, + parameters, } } } /// Pure component association parameters. #[derive(Serialize, Deserialize, Clone, Copy)] -pub struct AssociationRecord { - /// Association volume parameter - #[serde(skip_serializing_if = "f64::is_zero")] - #[serde(default)] - pub kappa_ab: f64, - /// Association energy parameter in units of Kelvin - #[serde(skip_serializing_if = "f64::is_zero")] - #[serde(default)] - pub epsilon_k_ab: f64, +pub struct AssociationRecord { + #[serde(flatten)] + pub parameters: A, /// \# of association sites of type A #[serde(skip_serializing_if = "f64::is_zero")] #[serde(default)] @@ -64,11 +52,10 @@ pub struct AssociationRecord { pub nc: f64, } -impl AssociationRecord { - pub fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { +impl AssociationRecord { + pub fn new(parameters: A, na: f64, nb: f64, nc: f64) -> Self { Self { - kappa_ab, - epsilon_k_ab, + parameters, na, nb, nc, @@ -76,10 +63,9 @@ impl AssociationRecord { } } -impl fmt::Display for AssociationRecord { +impl fmt::Display for AssociationRecord { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "AssociationRecord(kappa_ab={}", self.kappa_ab)?; - write!(f, ", epsilon_k_ab={}", self.epsilon_k_ab)?; + write!(f, "AssociationRecord(parameters={}", self.parameters)?; if self.na > 0.0 { write!(f, ", na={}", self.na)?; } @@ -95,13 +81,10 @@ impl fmt::Display for AssociationRecord { /// Binary association parameters. #[derive(Serialize, Deserialize, Clone, Copy)] -pub struct BinaryAssociationRecord { - /// Cross-association association volume parameter. - #[serde(skip_serializing_if = "Option::is_none")] - pub kappa_ab: Option, - /// Cross-association energy parameter. - #[serde(skip_serializing_if = "Option::is_none")] - pub epsilon_k_ab: Option, +pub struct BinaryAssociationRecord { + // Binary association parameters + #[serde(flatten)] + pub parameters: B, /// Indices of sites that the record refers to. #[serde(skip_serializing_if = "is_default_site_indices")] #[serde(default)] @@ -112,15 +95,10 @@ fn is_default_site_indices([i, j]: &[usize; 2]) -> bool { *i == 0 && *j == 0 } -impl BinaryAssociationRecord { - pub fn new( - kappa_ab: Option, - epsilon_k_ab: Option, - site_indices: Option<[usize; 2]>, - ) -> Self { +impl BinaryAssociationRecord { + pub fn new(parameters: B, site_indices: Option<[usize; 2]>) -> Self { Self { - kappa_ab, - epsilon_k_ab, + parameters, site_indices: site_indices.unwrap_or_default(), } } @@ -129,22 +107,19 @@ impl BinaryAssociationRecord { /// Parameter set required for the SAFT association Helmoltz energy /// contribution and functional. #[derive(Clone)] -pub struct AssociationParameters { +pub struct AssociationParameters { component_index: Array1, - sites_a: Array1, - sites_b: Array1, - sites_c: Array1, - pub sigma3_kappa_ab: Array2, - pub sigma3_kappa_cc: Array2, - pub epsilon_k_ab: Array2, - pub epsilon_k_cc: Array2, + sites_a: Array1>, + sites_b: Array1>, + sites_c: Array1>, + parameters_ab: Array2, + parameters_cc: Array2, } -impl AssociationParameters { +impl AssociationParameters

{ pub fn new( - records: &[Vec], - sigma: &Array1, - binary_records: &[((usize, usize), BinaryAssociationRecord)], + records: &[Vec>], + binary_records: &[([usize; 2], BinaryAssociationRecord)], component_index: Option<&Array1>, ) -> Self { let mut sites_a = Vec::new(); @@ -154,31 +129,13 @@ impl AssociationParameters { for (i, record) in records.iter().enumerate() { for (s, site) in record.iter().enumerate() { if site.na > 0.0 { - sites_a.push(AssociationSite::new( - i, - s, - site.na, - site.kappa_ab, - site.epsilon_k_ab, - )); + sites_a.push(AssociationSite::new(i, s, site.na, site.parameters)); } if site.nb > 0.0 { - sites_b.push(AssociationSite::new( - i, - s, - site.nb, - site.kappa_ab, - site.epsilon_k_ab, - )); + sites_b.push(AssociationSite::new(i, s, site.nb, site.parameters)); } if site.nc > 0.0 { - sites_c.push(AssociationSite::new( - i, - s, - site.nc, - site.kappa_ab, - site.epsilon_k_ab, - )); + sites_c.push(AssociationSite::new(i, s, site.nc, site.parameters)); } } } @@ -201,49 +158,24 @@ impl AssociationParameters { .map(|(i, site)| ((site.assoc_comp, site.site_index), i)) .collect(); - let mut sigma3_kappa_ab = - Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { - (sigma[sites_a[i].assoc_comp] * sigma[sites_b[j].assoc_comp]).powf(1.5) - * (sites_a[i].kappa_ab * sites_b[j].kappa_ab).sqrt() - }); - let mut sigma3_kappa_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { - (sigma[sites_c[i].assoc_comp] * sigma[sites_c[j].assoc_comp]).powf(1.5) - * (sites_c[i].kappa_ab * sites_c[j].kappa_ab).sqrt() + let mut parameters_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { + P::combining_rule(sites_a[i].parameters, sites_b[j].parameters) }); - let mut epsilon_k_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { - 0.5 * (sites_a[i].epsilon_k_ab + sites_b[j].epsilon_k_ab) - }); - let mut epsilon_k_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { - 0.5 * (sites_c[i].epsilon_k_ab + sites_c[j].epsilon_k_ab) + let mut parameters_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { + P::combining_rule(sites_c[i].parameters, sites_c[j].parameters) }); - for &((i, j), record) in binary_records.iter() { + for &([i, j], record) in binary_records.iter() { let [a, b] = record.site_indices; if let (Some(x), Some(y)) = (indices_a.get(&(i, a)), indices_b.get(&(j, b))) { - if let Some(epsilon_k_aibj) = record.epsilon_k_ab { - epsilon_k_ab[[*x, *y]] = epsilon_k_aibj; - } - if let Some(kappa_aibj) = record.kappa_ab { - sigma3_kappa_ab[[*x, *y]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj; - } + P::update_binary(&mut parameters_ab[[*x, *y]], record.parameters); } if let (Some(y), Some(x)) = (indices_b.get(&(i, a)), indices_a.get(&(j, b))) { - if let Some(epsilon_k_aibj) = record.epsilon_k_ab { - epsilon_k_ab[[*x, *y]] = epsilon_k_aibj; - } - if let Some(kappa_aibj) = record.kappa_ab { - sigma3_kappa_ab[[*x, *y]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj; - } + P::update_binary(&mut parameters_ab[[*x, *y]], record.parameters); } if let (Some(x), Some(y)) = (indices_c.get(&(i, a)), indices_c.get(&(j, b))) { - if let Some(epsilon_k_aibj) = record.epsilon_k_ab { - epsilon_k_cc[[*x, *y]] = epsilon_k_aibj; - epsilon_k_cc[[*y, *x]] = epsilon_k_aibj; - } - if let Some(kappa_aibj) = record.kappa_ab { - sigma3_kappa_cc[[*x, *y]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj; - sigma3_kappa_cc[[*y, *x]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj; - } + P::update_binary(&mut parameters_cc[[*x, *y]], record.parameters); + P::update_binary(&mut parameters_cc[[*y, *x]], record.parameters); } } @@ -254,10 +186,8 @@ impl AssociationParameters { sites_a: Array1::from_vec(sites_a), sites_b: Array1::from_vec(sites_b), sites_c: Array1::from_vec(sites_c), - sigma3_kappa_ab, - sigma3_kappa_cc, - epsilon_k_ab, - epsilon_k_cc, + parameters_ab, + parameters_cc, } } @@ -268,24 +198,24 @@ impl AssociationParameters { /// Implementation of the SAFT association Helmholtz energy /// contribution and functional. -pub struct Association

{ +pub struct Association { parameters: Arc

, - association_parameters: Arc, + association_parameters: Arc>, max_iter: usize, tol: f64, force_cross_association: bool, } -impl Association

{ +impl Association

{ pub fn new( parameters: &Arc

, - association_parameters: &AssociationParameters, + association_parameters: &Arc>, max_iter: usize, tol: f64, ) -> Self { Self { parameters: parameters.clone(), - association_parameters: Arc::new(association_parameters.clone()), + association_parameters: association_parameters.clone(), max_iter, tol, force_cross_association: false, @@ -294,7 +224,7 @@ impl Association

{ pub fn new_cross_association( parameters: &Arc

, - association_parameters: &AssociationParameters, + association_parameters: &Arc>, max_iter: usize, tol: f64, ) -> Self { @@ -302,48 +232,38 @@ impl Association

{ res.force_cross_association = true; res } +} + +pub trait AssociationStrength: HardSphereProperties { + type Record: Copy; + type BinaryRecord: Copy; fn association_strength + Copy>( &self, temperature: D, - diameter: &Array1, - n2: D, - n3i: D, - xi: D, - ) -> [Array2; 2] { - let p = &self.association_parameters; - let delta_ab = Array2::from_shape_fn([p.sites_a.len(), p.sites_b.len()], |(i, j)| { - let di = diameter[p.sites_a[i].assoc_comp]; - let dj = diameter[p.sites_b[j].assoc_comp]; - let k = di * dj / (di + dj) * (n2 * n3i); - n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) - * p.sigma3_kappa_ab[(i, j)] - * (temperature.recip() * p.epsilon_k_ab[(i, j)]).exp_m1() - }); - let delta_cc = Array2::from_shape_fn([p.sites_c.len(); 2], |(i, j)| { - let di = diameter[p.sites_c[i].assoc_comp]; - let dj = diameter[p.sites_c[j].assoc_comp]; - let k = di * dj / (di + dj) * (n2 * n3i); - n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) - * p.sigma3_kappa_cc[(i, j)] - * (temperature.recip() * p.epsilon_k_cc[(i, j)]).exp_m1() - }); - [delta_ab, delta_cc] - } + comp_i: usize, + comp_j: usize, + assoc_ij: Self::Record, + ) -> D; + + fn combining_rule(parameters_i: Self::Record, parameters_j: Self::Record) -> Self::Record; + + fn update_binary(_parameters_ij: &mut Self::Record, _binary_parameters: Self::BinaryRecord) {} } -impl Association

{ +impl Association

{ #[inline] pub fn helmholtz_energy + Copy>( &self, state: &StateHD, diameter: &Array1, ) -> D { - let p: &P = &self.parameters; let a = &self.association_parameters; // auxiliary variables - let [zeta2, n3] = p.zeta(state.temperature, &state.partial_density, [2, 3]); + let [zeta2, n3] = self + .parameters + .zeta(state.temperature, &state.partial_density, [2, 3]); let n2 = zeta2 * 6.0; let n3i = (-n3 + 1.0).recip(); @@ -387,15 +307,52 @@ impl Association

{ } } } + + fn association_strength + Copy>( + &self, + temperature: D, + diameter: &Array1, + n2: D, + n3i: D, + xi: D, + ) -> [Array2; 2] { + let p = &self.association_parameters; + + let delta_ab = Array2::from_shape_fn([p.sites_a.len(), p.sites_b.len()], |(i, j)| { + let di = diameter[p.sites_a[i].assoc_comp]; + let dj = diameter[p.sites_b[j].assoc_comp]; + let k = di * dj / (di + dj) * (n2 * n3i); + n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) + * self.parameters.association_strength( + temperature, + p.sites_a[i].assoc_comp, + p.sites_b[j].assoc_comp, + p.parameters_ab[(i, j)], + ) + }); + let delta_cc = Array2::from_shape_fn([p.sites_c.len(); 2], |(i, j)| { + let di = diameter[p.sites_c[i].assoc_comp]; + let dj = diameter[p.sites_c[j].assoc_comp]; + let k = di * dj / (di + dj) * (n2 * n3i); + n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) + * self.parameters.association_strength( + temperature, + p.sites_c[i].assoc_comp, + p.sites_c[j].assoc_comp, + p.parameters_cc[(i, j)], + ) + }); + [delta_ab, delta_cc] + } } -impl

fmt::Display for Association

{ +impl fmt::Display for Association

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

{ +impl Association

{ fn helmholtz_energy_ab_analytic + Copy>( &self, state: &StateHD, @@ -563,73 +520,77 @@ impl Association

{ } #[cfg(test)] -mod tests { +#[cfg(feature = "pcsaft")] +mod tests_pcsaft { use super::*; + use crate::hard_sphere::HardSphereProperties; + use crate::pcsaft::parameters::utils::water_parameters; + use crate::pcsaft::parameters::{PcSaftAssociationRecord, PcSaftBinaryAssociationRecord}; + use crate::pcsaft::PcSaftParameters; + use approx::assert_relative_eq; + use feos_core::parameter::{Parameter, ParameterError}; + + fn record( + kappa_ab: f64, + epsilon_k_ab: f64, + na: f64, + nb: f64, + ) -> AssociationRecord { + let pcsaft = PcSaftAssociationRecord::new(kappa_ab, epsilon_k_ab); + AssociationRecord::new(pcsaft, na, nb, 0.0) + } + + fn binary_record( + kappa_ab: f64, + epsilon_k_ab: f64, + indices: Option<[usize; 2]>, + ) -> BinaryAssociationRecord { + let pcsaft = PcSaftBinaryAssociationRecord::new(Some(kappa_ab), Some(epsilon_k_ab)); + BinaryAssociationRecord::new(pcsaft, indices) + } #[test] fn test_binary_parameters() { - let comp1 = vec![AssociationRecord::new(0.1, 2500., 1.0, 1.0, 0.0)]; - let comp2 = vec![AssociationRecord::new(0.2, 1500., 1.0, 1.0, 0.0)]; - let comp3 = vec![AssociationRecord::new(0.3, 500., 0.0, 1.0, 0.0)]; - let comp4 = vec![ - AssociationRecord::new(0.3, 1000., 1.0, 0.0, 0.0), - AssociationRecord::new(0.3, 2000., 0.0, 1.0, 0.0), - ]; + let comp1 = vec![record(0.1, 2500., 1.0, 1.0)]; + let comp2 = vec![record(0.2, 1500., 1.0, 1.0)]; + let comp3 = vec![record(0.3, 500., 0.0, 1.0)]; + let comp4 = vec![record(0.3, 1000., 1.0, 0.0), record(0.3, 2000., 0.0, 1.0)]; let records = [comp1, comp2, comp3, comp4]; - let sigma = arr1(&[3.0, 3.0, 3.0, 3.0]); let binary = [ - ( - (0, 1), - BinaryAssociationRecord::new(Some(3.5), Some(1234.), Some([0, 0])), - ), - ( - (0, 2), - BinaryAssociationRecord::new(Some(3.5), Some(3140.), Some([0, 0])), - ), - ( - (1, 3), - BinaryAssociationRecord::new(Some(3.5), Some(3333.), Some([0, 1])), - ), + ([0, 1], binary_record(3.5, 1234., Some([0, 0]))), + ([0, 2], binary_record(3.5, 3140., Some([0, 0]))), + ([1, 3], binary_record(3.5, 3333., Some([0, 1]))), ]; - let assoc = AssociationParameters::new(&records, &sigma, &binary, None); - println!("{}", assoc.epsilon_k_ab); + let assoc = AssociationParameters::::new(&records, &binary, None); + println!("{}", assoc.parameters_ab); let epsilon_k_ab = arr2(&[ [2500., 1234., 3140., 2250.], [1234., 1500., 1000., 3333.], [1750., 1250., 750., 1500.], ]); - assert_eq!(assoc.epsilon_k_ab, epsilon_k_ab); + assert_eq!(assoc.parameters_ab.mapv(|p| p.epsilon_k_ab), epsilon_k_ab); } #[test] fn test_induced_association() { - let comp1 = vec![AssociationRecord::new(0.1, 2500., 1.0, 1.0, 0.0)]; - let comp2 = vec![AssociationRecord::new(0.1, -500., 0.0, 1.0, 0.0)]; - let comp3 = vec![AssociationRecord::new(0.0, 0.0, 0.0, 1.0, 0.0)]; - let sigma = arr1(&[3.0, 3.5]); - let binary = [( - (0, 1), - BinaryAssociationRecord::new(Some(0.1), Some(1000.), None), - )]; - let assoc1 = AssociationParameters::new(&[comp1.clone(), comp2], &sigma, &[], None); - let assoc2 = AssociationParameters::new(&[comp1, comp3], &sigma, &binary, None); - println!("{}", assoc1.epsilon_k_ab); - println!("{}", assoc2.epsilon_k_ab); - assert_eq!(assoc1.epsilon_k_ab, assoc2.epsilon_k_ab); - println!("{}", assoc1.sigma3_kappa_ab); - println!("{}", assoc2.sigma3_kappa_ab); - assert_eq!(assoc1.sigma3_kappa_ab, assoc2.sigma3_kappa_ab); + let comp1 = vec![record(0.1, 2500., 1.0, 1.0)]; + let comp2 = vec![record(0.1, -500., 0.0, 1.0)]; + let comp3 = vec![record(0.0, 0.0, 0.0, 1.0)]; + let binary = [([0, 1], binary_record(0.1, 1000., None))]; + let assoc1 = + AssociationParameters::::new(&[comp1.clone(), comp2], &[], None); + let assoc2 = AssociationParameters::::new(&[comp1, comp3], &binary, None); + println!("{}", assoc1.parameters_ab); + println!("{}", assoc2.parameters_ab); + assert_eq!( + assoc1.parameters_ab.mapv(|p| p.epsilon_k_ab), + assoc2.parameters_ab.mapv(|p| p.epsilon_k_ab) + ); + assert_eq!( + assoc1.parameters_ab.mapv(|p| p.kappa_ab), + assoc2.parameters_ab.mapv(|p| p.kappa_ab) + ); } -} - -#[cfg(test)] -#[cfg(feature = "pcsaft")] -mod tests_pcsaft { - use super::*; - use crate::pcsaft::parameters::utils::water_parameters; - use crate::pcsaft::PcSaftParameters; - use approx::assert_relative_eq; - use feos_core::parameter::{Parameter, ParameterError}; #[test] fn helmholtz_energy() { diff --git a/src/association/python.rs b/src/association/python.rs deleted file mode 100644 index 9744f5552..000000000 --- a/src/association/python.rs +++ /dev/null @@ -1,49 +0,0 @@ -use super::AssociationRecord; -use feos_core::impl_json_handling; -use feos_core::parameter::ParameterError; -use pyo3::prelude::*; - -/// Pure component association parameters -#[pyclass(name = "AssociationRecord")] -#[derive(Clone)] -pub struct PyAssociationRecord(pub AssociationRecord); - -#[pymethods] -impl PyAssociationRecord { - #[new] - #[pyo3(signature = (kappa_ab, epsilon_k_ab, na=0.0, nb=0.0, nc=0.0))] - fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { - Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc)) - } - - #[getter] - fn get_kappa_ab(&self) -> f64 { - self.0.kappa_ab - } - - #[getter] - fn get_epsilon_k_ab(&self) -> f64 { - self.0.epsilon_k_ab - } - - #[getter] - fn get_na(&self) -> f64 { - self.0.na - } - - #[getter] - fn get_nb(&self) -> f64 { - self.0.nb - } - - #[getter] - fn get_nc(&self) -> f64 { - self.0.nc - } - - fn __repr__(&self) -> PyResult { - Ok(self.0.to_string()) - } -} - -impl_json_handling!(PyAssociationRecord); diff --git a/src/epcsaft/mod.rs b/src/epcsaft/mod.rs index 834f53c22..5a0d05662 100644 --- a/src/epcsaft/mod.rs +++ b/src/epcsaft/mod.rs @@ -1,8 +1,4 @@ -//! Electrolyte Perturbed-Chain Statistical Associating Fluid Theory (e12PC-SAFT) -//! - -#![warn(clippy::all)] -#![allow(clippy::too_many_arguments)] +//! Electrolyte Perturbed-Chain Statistical Associating Fluid Theory (ePC-SAFT) mod eos; pub(crate) mod parameters; diff --git a/src/epcsaft/parameters.rs b/src/epcsaft/parameters.rs index 596e6e923..21ef0407c 100644 --- a/src/epcsaft/parameters.rs +++ b/src/epcsaft/parameters.rs @@ -1,13 +1,15 @@ -use crate::association::{AssociationParameters, AssociationRecord, BinaryAssociationRecord}; +use crate::association::{ + AssociationParameters, AssociationRecord, AssociationStrength, BinaryAssociationRecord, +}; use crate::hard_sphere::{HardSphereProperties, MonomerShape}; use feos_core::parameter::{FromSegments, Parameter, ParameterError, PureRecord}; -use feos_core::si::{JOULE, KB, KELVIN}; use ndarray::{Array, Array1, Array2}; use num_dual::DualNum; use num_traits::Zero; use serde::{Deserialize, Serialize}; use std::collections::HashMap; use std::fmt::Write; +use std::sync::Arc; use crate::epcsaft::eos::permittivity::PermittivityRecord; @@ -20,16 +22,10 @@ pub struct ElectrolytePcSaftRecord { pub sigma: f64, /// Energetic parameter in units of Kelvin pub epsilon_k: f64, - /// Dipole moment in units of Debye - #[serde(skip_serializing_if = "Option::is_none")] - pub mu: Option, - /// Quadrupole moment in units of Debye - #[serde(skip_serializing_if = "Option::is_none")] - pub q: Option, /// Association parameters #[serde(flatten)] #[serde(skip_serializing_if = "Option::is_none")] - pub association_record: Option, + pub association_record: Option>, #[serde(default)] #[serde(skip_serializing_if = "Option::is_none")] pub z: Option, @@ -51,21 +47,13 @@ impl FromSegments for ElectrolytePcSaftRecord { z += s.z.unwrap_or(0.0); }); - let q = segments - .iter() - .filter_map(|(s, n)| s.q.map(|q| q * n)) - .reduce(|a, b| a + b); - let mu = segments - .iter() - .filter_map(|(s, n)| s.mu.map(|mu| mu * n)) - .reduce(|a, b| a + b); let association_record = segments .iter() .filter_map(|(s, n)| { s.association_record.as_ref().map(|record| { [ - record.kappa_ab * n, - record.epsilon_k_ab * n, + record.parameters.kappa_ab * n, + record.parameters.epsilon_k_ab * n, record.na * n, record.nb * n, record.nc * n, @@ -82,15 +70,18 @@ impl FromSegments for ElectrolytePcSaftRecord { ] }) .map(|[kappa_ab, epsilon_k_ab, na, nb, nc]| { - AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc) + AssociationRecord::new( + ElectrolytePcSaftAssociationRecord::new(kappa_ab, epsilon_k_ab), + na, + nb, + nc, + ) }); Ok(Self { m, sigma: (sigma3 / m).cbrt(), epsilon_k: epsilon_k / m, - mu, - q, association_record, z: Some(z), permittivity_record: None, @@ -115,12 +106,6 @@ impl std::fmt::Display for ElectrolytePcSaftRecord { write!(f, "ElectrolytePcSaftRecord(m={}", self.m)?; write!(f, ", sigma={}", self.sigma)?; write!(f, ", epsilon_k={}", self.epsilon_k)?; - if let Some(n) = &self.mu { - write!(f, ", mu={}", n)?; - } - if let Some(n) = &self.q { - write!(f, ", q={}", n)?; - } if let Some(n) = &self.association_record { write!(f, ", association_record={}", n)?; } @@ -139,8 +124,6 @@ impl ElectrolytePcSaftRecord { m: f64, sigma: f64, epsilon_k: f64, - mu: Option, - q: Option, kappa_ab: Option, epsilon_k_ab: Option, na: Option, @@ -149,28 +132,21 @@ impl ElectrolytePcSaftRecord { z: Option, permittivity_record: Option, ) -> ElectrolytePcSaftRecord { - let association_record = if kappa_ab.is_none() - && epsilon_k_ab.is_none() - && na.is_none() - && nb.is_none() - && nc.is_none() - { - None - } else { - Some(AssociationRecord::new( - kappa_ab.unwrap_or_default(), - epsilon_k_ab.unwrap_or_default(), - na.unwrap_or_default(), - nb.unwrap_or_default(), - nc.unwrap_or_default(), - )) - }; + let association_record = + if let (Some(kappa_ab), Some(epsilon_k_ab)) = (kappa_ab, epsilon_k_ab) { + Some(AssociationRecord::new( + ElectrolytePcSaftAssociationRecord::new(kappa_ab, epsilon_k_ab), + na.unwrap_or_default(), + nb.unwrap_or_default(), + nc.unwrap_or_default(), + )) + } else { + None + }; ElectrolytePcSaftRecord { m, sigma, epsilon_k, - mu, - q, association_record, z, permittivity_record, @@ -178,6 +154,34 @@ impl ElectrolytePcSaftRecord { } } +#[derive(Serialize, Deserialize, Clone, Copy, Default)] +pub struct ElectrolytePcSaftAssociationRecord { + /// Association volume parameter + pub kappa_ab: f64, + /// Association energy parameter in units of Kelvin + pub epsilon_k_ab: f64, +} + +impl ElectrolytePcSaftAssociationRecord { + pub fn new(kappa_ab: f64, epsilon_k_ab: f64) -> Self { + Self { + kappa_ab, + epsilon_k_ab, + } + } +} + +impl std::fmt::Display for ElectrolytePcSaftAssociationRecord { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!( + f, + "ElectrolytePcSaftAssociationRecord(kappa_ab={}", + self.kappa_ab + )?; + write!(f, ", epsilon_k_ab={})", self.epsilon_k_ab) + } +} + #[derive(Serialize, Deserialize, Clone, Default)] pub struct ElectrolytePcSaftBinaryRecord { /// Binary dispersion interaction parameter @@ -185,7 +189,7 @@ pub struct ElectrolytePcSaftBinaryRecord { pub k_ij: Vec, /// Binary association parameters #[serde(flatten)] - association: Option, + association: Option>, } impl ElectrolytePcSaftBinaryRecord { @@ -194,7 +198,10 @@ impl ElectrolytePcSaftBinaryRecord { let association = if kappa_ab.is_none() && epsilon_k_ab.is_none() { None } else { - Some(BinaryAssociationRecord::new(kappa_ab, epsilon_k_ab, None)) + Some(BinaryAssociationRecord::new( + ElectrolytePcSaftBinaryAssociationRecord::new(kappa_ab, epsilon_k_ab), + None, + )) }; Self { k_ij, association } } @@ -243,10 +250,10 @@ impl std::fmt::Display for ElectrolytePcSaftBinaryRecord { tokens.push(format!(", k_ij_3={})", self.k_ij[3])); } if let Some(association) = self.association { - if let Some(kappa_ab) = association.kappa_ab { + if let Some(kappa_ab) = association.parameters.kappa_ab { tokens.push(format!("kappa_ab={}", kappa_ab)); } - if let Some(epsilon_k_ab) = association.epsilon_k_ab { + if let Some(epsilon_k_ab) = association.parameters.epsilon_k_ab { tokens.push(format!("epsilon_k_ab={}", epsilon_k_ab)); } } @@ -254,27 +261,38 @@ impl std::fmt::Display for ElectrolytePcSaftBinaryRecord { } } +#[derive(Serialize, Deserialize, Clone, Copy, Default)] +pub struct ElectrolytePcSaftBinaryAssociationRecord { + /// Cross-association association volume parameter. + #[serde(skip_serializing_if = "Option::is_none")] + pub kappa_ab: Option, + /// Cross-association energy parameter. + #[serde(skip_serializing_if = "Option::is_none")] + pub epsilon_k_ab: Option, +} + +impl ElectrolytePcSaftBinaryAssociationRecord { + pub fn new(kappa_ab: Option, epsilon_k_ab: Option) -> Self { + Self { + kappa_ab, + epsilon_k_ab, + } + } +} + pub struct ElectrolytePcSaftParameters { pub molarweight: Array1, pub m: Array1, pub sigma: Array1, pub epsilon_k: Array1, - pub mu: Array1, - pub q: Array1, - pub mu2: Array1, - pub q2: Array1, - pub association: AssociationParameters, + pub association: Arc>, pub z: Array1, pub k_ij: Array2>, pub sigma_ij: Array2, pub e_k_ij: Array2, - pub ndipole: usize, - pub nquadpole: usize, pub nionic: usize, pub nsolvent: usize, pub water_sigma_t_comp: Option, - pub dipole_comp: Array1, - pub quadpole_comp: Array1, pub ionic_comp: Array1, pub solvent_comp: Array1, pub permittivity: Array1>, @@ -323,8 +341,6 @@ impl Parameter for ElectrolytePcSaftParameters { let mut m = Array::zeros(n); let mut sigma = Array::zeros(n); let mut epsilon_k = Array::zeros(n); - let mut mu = Array::zeros(n); - let mut q = Array::zeros(n); let mut z = Array::zeros(n); let mut association_records = Vec::with_capacity(n); let mut water_sigma_t_comp = None; @@ -337,16 +353,14 @@ impl Parameter for ElectrolytePcSaftParameters { m[i] = r.m; sigma[i] = r.sigma; epsilon_k[i] = r.epsilon_k; - mu[i] = r.mu.unwrap_or(0.0); - q[i] = r.q.unwrap_or(0.0); z[i] = r.z.unwrap_or(0.0); association_records.push(r.association_record.into_iter().collect()); molarweight[i] = record.molarweight; // check if component i is water with temperature-dependent sigma if (m[i] * 1000.0).round() / 1000.0 == 1.205 && epsilon_k[i].round() == 354.0 { if let Some(record) = r.association_record { - if (record.kappa_ab * 1000.0).round() / 1000.0 == 0.045 - && record.epsilon_k_ab.round() == 2426.0 + if (record.parameters.kappa_ab * 1000.0).round() / 1000.0 == 0.045 + && record.parameters.epsilon_k_ab.round() == 2426.0 { water_sigma_t_comp = Some(i); } @@ -354,34 +368,15 @@ impl Parameter for ElectrolytePcSaftParameters { } } - let mu2 = &mu * &mu / (&m * &sigma * &sigma * &sigma * &epsilon_k) - * 1e-19 - * (JOULE / KELVIN / KB).into_value(); - let q2 = &q * &q / (&m * &sigma.mapv(|s| s.powi(5)) * &epsilon_k) - * 1e-19 - * (JOULE / KELVIN / KB).into_value(); - let dipole_comp: Array1 = mu2 - .iter() - .enumerate() - .filter_map(|(i, &mu2)| (mu2.abs() > 0.0).then_some(i)) - .collect(); - let ndipole = dipole_comp.len(); - let quadpole_comp: Array1 = q2 - .iter() - .enumerate() - .filter_map(|(i, &q2)| (q2.abs() > 0.0).then_some(i)) - .collect(); - let nquadpole = quadpole_comp.len(); - let binary_association: Vec<_> = binary_records .iter() .flat_map(|r| { r.indexed_iter() - .filter_map(|(i, record)| record.association.map(|r| (i, r))) + .filter_map(|((i, j), record)| record.association.map(|r| ([i, j], r))) }) .collect(); let association = - AssociationParameters::new(&association_records, &sigma, &binary_association, None); + AssociationParameters::new(&association_records, &binary_association, None); let ionic_comp: Array1 = z .iter() @@ -504,21 +499,13 @@ impl Parameter for ElectrolytePcSaftParameters { m, sigma, epsilon_k, - mu, - q, - mu2, - q2, - association, + association: Arc::new(association), z, k_ij, sigma_ij, e_k_ij, - ndipole, - nquadpole, nionic, nsolvent, - dipole_comp, - quadpole_comp, ionic_comp, solvent_comp, water_sigma_t_comp, @@ -558,35 +545,73 @@ impl HardSphereProperties for ElectrolytePcSaftParameters { } } +impl AssociationStrength for ElectrolytePcSaftParameters { + type Record = ElectrolytePcSaftAssociationRecord; + type BinaryRecord = ElectrolytePcSaftBinaryAssociationRecord; + + fn association_strength + Copy>( + &self, + temperature: D, + comp_i: usize, + comp_j: usize, + assoc_ij: Self::Record, + ) -> D { + let sigma_t = self.sigma_t(temperature); + let si = sigma_t[comp_i]; + let sj = sigma_t[comp_j]; + (temperature.recip() * assoc_ij.epsilon_k_ab).exp_m1() + * assoc_ij.kappa_ab + * (si * sj).powf(1.5) + } + + fn combining_rule(parameters_i: Self::Record, parameters_j: Self::Record) -> Self::Record { + Self::Record { + kappa_ab: (parameters_i.kappa_ab * parameters_j.kappa_ab).sqrt(), + epsilon_k_ab: 0.5 * (parameters_i.epsilon_k_ab + parameters_j.epsilon_k_ab), + } + } + + fn update_binary(parameters_ij: &mut Self::Record, binary_parameters: Self::BinaryRecord) { + if let Some(kappa_ab) = binary_parameters.kappa_ab { + parameters_ij.kappa_ab = kappa_ab + } + if let Some(epsilon_k_ab) = binary_parameters.epsilon_k_ab { + parameters_ij.epsilon_k_ab = epsilon_k_ab + } + } +} + impl ElectrolytePcSaftParameters { pub fn to_markdown(&self) -> String { let mut output = String::new(); let o = &mut output; write!( o, - "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$\\mu$|$Q$|$z$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|-|" + "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$z$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); for (i, record) in self.pure_records.iter().enumerate() { let component = record.identifier.name.clone(); let component = component.unwrap_or(format!("Component {}", i + 1)); - let association = record - .model_record - .association_record - .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, 0.0, 0.0, 0.0)); + let association = record.model_record.association_record.unwrap_or_else(|| { + AssociationRecord::new( + ElectrolytePcSaftAssociationRecord::new(0.0, 0.0), + 0.0, + 0.0, + 0.0, + ) + }); write!( o, - "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", + "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", component, record.molarweight, record.model_record.m, record.model_record.sigma, record.model_record.epsilon_k, - record.model_record.mu.unwrap_or(0.0), - record.model_record.q.unwrap_or(0.0), record.model_record.z.unwrap_or(0.0), - association.kappa_ab, - association.epsilon_k_ab, + association.parameters.kappa_ab, + association.parameters.epsilon_k_ab, association.na, association.nb, association.nc @@ -598,7 +623,6 @@ impl ElectrolytePcSaftParameters { } } -#[allow(dead_code)] #[cfg(test)] pub mod utils { use feos_core::parameter::{BinaryRecord, Identifier}; @@ -629,30 +653,6 @@ pub mod utils { Arc::new(ElectrolytePcSaftParameters::new_pure(propane_record).unwrap()) } - pub fn carbon_dioxide_parameters() -> Arc { - let co2_json = r#" - { - "identifier": { - "cas": "124-38-9", - "name": "carbon-dioxide", - "iupac_name": "carbon dioxide", - "smiles": "O=C=O", - "inchi": "InChI=1/CO2/c2-1-3", - "formula": "CO2" - }, - "molarweight": 44.0098, - "model_record": { - "m": 1.5131, - "sigma": 3.1869, - "epsilon_k": 163.333, - "q": 4.4 - } - }"#; - let co2_record: PureRecord = - serde_json::from_str(co2_json).expect("Unable to parse json."); - Arc::new(ElectrolytePcSaftParameters::new_pure(co2_record).unwrap()) - } - pub fn butane_parameters() -> Arc { let butane_json = r#" { @@ -676,55 +676,6 @@ pub mod utils { Arc::new(ElectrolytePcSaftParameters::new_pure(butane_record).unwrap()) } - pub fn dme_parameters() -> Arc { - let dme_json = r#" - { - "identifier": { - "cas": "115-10-6", - "name": "dimethyl-ether", - "iupac_name": "methoxymethane", - "smiles": "COC", - "inchi": "InChI=1/C2H6O/c1-3-2/h1-2H3", - "formula": "C2H6O" - }, - "model_record": { - "m": 2.2634, - "sigma": 3.2723, - "epsilon_k": 210.29, - "mu": 1.3 - }, - "molarweight": 46.0688 - }"#; - let dme_record: PureRecord = - serde_json::from_str(dme_json).expect("Unable to parse json."); - Arc::new(ElectrolytePcSaftParameters::new_pure(dme_record).unwrap()) - } - - pub fn water_parameters_sigma_t() -> Arc { - let water_json = r#" - { - "identifier": { - "cas": "7732-18-5", - "name": "water_np_sigma_t", - "iupac_name": "oxidane", - "smiles": "O", - "inchi": "InChI=1/H2O/h1H2", - "formula": "H2O" - }, - "model_record": { - "m": 1.2047, - "sigma": 2.7927, - "epsilon_k": 353.95, - "kappa_ab": 0.04509, - "epsilon_k_ab": 2425.7 - }, - "molarweight": 18.0152 - }"#; - let water_record: PureRecord = - serde_json::from_str(water_json).expect("Unable to parse json."); - Arc::new(ElectrolytePcSaftParameters::new_pure(water_record).unwrap()) - } - pub fn water_nacl_parameters_perturb() -> Arc { // Water parameters from Held et al. (2014), originally from Fuchs et al. (2006) let pure_json = r#"[ @@ -1092,48 +1043,6 @@ pub mod utils { Arc::new(ElectrolytePcSaftParameters::new_pure(water_record).unwrap()) } - pub fn dme_co2_parameters() -> Arc { - let binary_json = r#"[ - { - "identifier": { - "cas": "115-10-6", - "name": "dimethyl-ether", - "iupac_name": "methoxymethane", - "smiles": "COC", - "inchi": "InChI=1/C2H6O/c1-3-2/h1-2H3", - "formula": "C2H6O" - }, - "molarweight": 46.0688, - "model_record": { - "m": 2.2634, - "sigma": 3.2723, - "epsilon_k": 210.29, - "mu": 1.3 - } - }, - { - "identifier": { - "cas": "124-38-9", - "name": "carbon-dioxide", - "iupac_name": "carbon dioxide", - "smiles": "O=C=O", - "inchi": "InChI=1/CO2/c2-1-3", - "formula": "CO2" - }, - "molarweight": 44.0098, - "model_record": { - "m": 1.5131, - "sigma": 3.1869, - "epsilon_k": 163.333, - "q": 4.4 - } - } - ]"#; - let binary_record: Vec> = - serde_json::from_str(binary_json).expect("Unable to parse json."); - Arc::new(ElectrolytePcSaftParameters::new_binary(binary_record, None).unwrap()) - } - pub fn propane_butane_parameters() -> Arc { let binary_json = r#"[ { diff --git a/src/epcsaft/python.rs b/src/epcsaft/python.rs index 9a9a0a415..e90fa92cb 100644 --- a/src/epcsaft/python.rs +++ b/src/epcsaft/python.rs @@ -25,10 +25,6 @@ use std::sync::Arc; /// Segment diameter in units of Angstrom. /// epsilon_k : float /// Energetic parameter in units of Kelvin. -/// mu : float, optional -/// Dipole moment in units of Debye. -/// q : float, optional -/// Quadrupole moment in units of Debye * Angstrom. /// kappa_ab : float, optional /// Association volume parameter. /// epsilon_k_ab : float, optional @@ -43,7 +39,6 @@ use std::sync::Arc; /// Charge of the electrolyte. /// permittivity_record : PyPermittivityRecord, optional /// Permittivity record. Defaults to `None`. - #[pyclass(name = "ElectrolytePcSaftRecord")] #[derive(Clone)] pub struct PyElectrolytePcSaftRecord(ElectrolytePcSaftRecord); @@ -52,14 +47,12 @@ pub struct PyElectrolytePcSaftRecord(ElectrolytePcSaftRecord); impl PyElectrolytePcSaftRecord { #[new] #[pyo3( - text_signature = "(m, sigma, epsilon_k, mu=None, q=None, kappa_ab=None, epsilon_k_ab=None, na=None, nb=None, nc=None, permittivity_record=None)" + text_signature = "(m, sigma, epsilon_k, kappa_ab=None, epsilon_k_ab=None, na=None, nb=None, nc=None, permittivity_record=None)" )] fn new( m: f64, sigma: f64, epsilon_k: f64, - mu: Option, - q: Option, kappa_ab: Option, epsilon_k_ab: Option, na: Option, @@ -68,23 +61,17 @@ impl PyElectrolytePcSaftRecord { z: Option, permittivity_record: Option, ) -> Self { - let perm = match permittivity_record { - Some(p) => Some(p.0), - None => None, - }; Self(ElectrolytePcSaftRecord::new( m, sigma, epsilon_k, - mu, - q, kappa_ab, epsilon_k_ab, na, nb, nc, z, - perm, + permittivity_record.map(|p| p.0), )) } @@ -105,12 +92,12 @@ impl PyElectrolytePcSaftRecord { #[getter] fn get_kappa_ab(&self) -> Option { - self.0.association_record.map(|a| a.kappa_ab) + self.0.association_record.map(|a| a.parameters.kappa_ab) } #[getter] fn get_epsilon_k_ab(&self) -> Option { - self.0.association_record.map(|a| a.epsilon_k_ab) + self.0.association_record.map(|a| a.parameters.epsilon_k_ab) } #[getter] @@ -212,7 +199,6 @@ impl PyPermittivityRecord { /// PermittivityRecord /// #[staticmethod] - #[allow(non_snake_case)] #[pyo3(text_signature = "(interpolation_points)")] pub fn from_experimental_data(interpolation_points: Vec<(f64, f64)>) -> Self { Self(PermittivityRecord::ExperimentalData { @@ -233,7 +219,6 @@ impl PyPermittivityRecord { /// PermittivityRecord /// #[staticmethod] - #[allow(non_snake_case)] #[pyo3( text_signature = "(dipole_scaling, polarizability_scaling, correlation_integral_parameter)" )] diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 65809f69f..4961906b0 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -1,5 +1,6 @@ use super::eos::GcPcSaftOptions; -use crate::association::Association; +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::si::{MolarWeight, GRAM, MOL}; @@ -143,6 +144,32 @@ impl HardSphereProperties for GcPcSaftFunctionalParameters { } } +impl AssociationStrength for GcPcSaftFunctionalParameters { + type Record = GcPcSaftAssociationRecord; + type BinaryRecord = (); + + fn association_strength + Copy>( + &self, + temperature: D, + comp_i: usize, + comp_j: usize, + assoc_ij: Self::Record, + ) -> D { + let si = self.sigma[comp_i]; + let sj = self.sigma[comp_j]; + (temperature.recip() * assoc_ij.epsilon_k_ab).exp_m1() + * assoc_ij.kappa_ab + * (si * sj).powf(1.5) + } + + fn combining_rule(parameters_i: Self::Record, parameters_j: Self::Record) -> Self::Record { + Self::Record { + kappa_ab: (parameters_i.kappa_ab * parameters_j.kappa_ab).sqrt(), + epsilon_k_ab: 0.5 * (parameters_i.epsilon_k_ab + parameters_j.epsilon_k_ab), + } + } +} + impl FluidParameters for GcPcSaftFunctional { fn epsilon_k_ff(&self) -> Array1 { self.parameters.epsilon_k.clone() diff --git a/src/gc_pcsaft/dft/parameter.rs b/src/gc_pcsaft/dft/parameter.rs index 68d644b0d..808a221af 100644 --- a/src/gc_pcsaft/dft/parameter.rs +++ b/src/gc_pcsaft/dft/parameter.rs @@ -7,6 +7,7 @@ use indexmap::IndexMap; use ndarray::{Array1, Array2}; use petgraph::dot::{Config, Dot}; use petgraph::graph::{Graph, UnGraph}; +use std::sync::Arc; /// psi Parameter for heterosegmented DFT (Mairhofer2018) const PSI_GC_DFT: f64 = 1.5357; @@ -20,7 +21,7 @@ pub struct GcPcSaftFunctionalParameters { pub sigma: Array1, pub epsilon_k: Array1, pub bonds: UnGraph<(), ()>, - pub association: AssociationParameters, + pub association: Arc>, pub psi_dft: Array1, pub k_ij: Array2, pub sigma_ij: Array2, @@ -128,7 +129,7 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { let sigma = Array1::from_vec(sigma); let component_index = Array1::from_vec(component_index); let association = - AssociationParameters::new(&association_records, &sigma, &[], Some(&component_index)); + AssociationParameters::new(&association_records, &[], Some(&component_index)); Ok(Self { molarweight, @@ -138,7 +139,7 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { sigma, epsilon_k: Array1::from_vec(epsilon_k), bonds, - association, + association: Arc::new(association), psi_dft: Array1::from_vec(psi_dft), k_ij, sigma_ij, diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs index ba365661b..5ba0311b7 100644 --- a/src/gc_pcsaft/eos/parameter.rs +++ b/src/gc_pcsaft/eos/parameter.rs @@ -1,5 +1,5 @@ -use crate::association::AssociationParameters; -use crate::gc_pcsaft::record::GcPcSaftRecord; +use crate::association::{AssociationParameters, AssociationStrength}; +use crate::gc_pcsaft::record::{GcPcSaftAssociationRecord, GcPcSaftRecord}; use crate::hard_sphere::{HardSphereProperties, MonomerShape}; use feos_core::parameter::{ BinaryRecord, ChemicalRecord, Identifier, ParameterError, ParameterHetero, SegmentCount, @@ -12,6 +12,7 @@ use num_dual::DualNum; use std::borrow::Cow; use std::collections::HashMap; use std::fmt::Write; +use std::sync::Arc; /// [ChemicalRecord] that is used as input for the gc-PC-SAFT equation of state. #[derive(Clone)] @@ -74,7 +75,7 @@ pub struct GcPcSaftEosParameters { pub epsilon_k: Array1, pub bonds: IndexMap<[usize; 2], f64>, - pub association: AssociationParameters, + pub association: Arc>, pub dipole_comp: Array1, mu: Array1, @@ -224,7 +225,7 @@ impl ParameterHetero for GcPcSaftEosParameters { let sigma = Array1::from_vec(sigma); let component_index = Array1::from_vec(component_index); let association = - AssociationParameters::new(&association_records, &sigma, &[], Some(&component_index)); + AssociationParameters::new(&association_records, &[], Some(&component_index)); Ok(Self { molarweight, @@ -235,7 +236,7 @@ impl ParameterHetero for GcPcSaftEosParameters { sigma, epsilon_k: Array1::from_vec(epsilon_k), bonds, - association, + association: Arc::new(association), dipole_comp: Array1::from_vec(dipole_comp), mu: Array1::from_vec(mu), mu2: Array1::from_vec(mu2), @@ -288,6 +289,32 @@ impl HardSphereProperties for GcPcSaftEosParameters { } } +impl AssociationStrength for GcPcSaftEosParameters { + type Record = GcPcSaftAssociationRecord; + type BinaryRecord = (); + + fn association_strength + Copy>( + &self, + temperature: D, + comp_i: usize, + comp_j: usize, + assoc_ij: Self::Record, + ) -> D { + let si = self.sigma[comp_i]; + let sj = self.sigma[comp_j]; + (temperature.recip() * assoc_ij.epsilon_k_ab).exp_m1() + * assoc_ij.kappa_ab + * (si * sj).powf(1.5) + } + + fn combining_rule(parameters_i: Self::Record, parameters_j: Self::Record) -> Self::Record { + Self::Record { + kappa_ab: (parameters_i.kappa_ab * parameters_j.kappa_ab).sqrt(), + epsilon_k_ab: 0.5 * (parameters_i.epsilon_k_ab + parameters_j.epsilon_k_ab), + } + } +} + impl GcPcSaftEosParameters { pub fn to_markdown(&self) -> String { let gorup_dict: HashMap<&String, &GcPcSaftRecord> = self @@ -329,7 +356,7 @@ impl GcPcSaftEosParameters { let association = if let Some(a) = record.association_record { format!( "{}|{}|{}|{}|{}", - a.kappa_ab, a.epsilon_k_ab, a.na, a.nb, a.nc + a.parameters.kappa_ab, a.parameters.epsilon_k_ab, a.na, a.nb, a.nc ) } else { "||||".to_string() @@ -408,14 +435,15 @@ impl std::fmt::Display for GcPcSaftEosParameters { #[cfg(test)] pub mod test { use super::*; - use crate::association::AssociationRecord; use feos_core::parameter::{ChemicalRecord, Identifier}; fn ch3() -> SegmentRecord { SegmentRecord::new( "CH3".into(), 15.0, - GcPcSaftRecord::new(0.77247, 3.6937, 181.49, None, None, None), + GcPcSaftRecord::new( + 0.77247, 3.6937, 181.49, None, None, None, None, None, None, None, + ), ) } @@ -423,7 +451,9 @@ pub mod test { SegmentRecord::new( "CH2".into(), 14.0, - GcPcSaftRecord::new(0.7912, 3.0207, 157.23, None, None, None), + GcPcSaftRecord::new( + 0.7912, 3.0207, 157.23, None, None, None, None, None, None, None, + ), ) } @@ -436,7 +466,11 @@ pub mod test { 2.7702, 334.29, None, - Some(AssociationRecord::new(0.009583, 2575.9, 1.0, 1.0, 0.0)), + Some(0.009583), + Some(2575.9), + Some(1.0), + Some(1.0), + None, None, ), ) diff --git a/src/gc_pcsaft/mod.rs b/src/gc_pcsaft/mod.rs index 09345e552..6f759ceb7 100644 --- a/src/gc_pcsaft/mod.rs +++ b/src/gc_pcsaft/mod.rs @@ -2,8 +2,6 @@ //! //! - [Gross et al. (2003)](https://doi.org/10.1021/ie020509y) //! - [Sauer et al. (2014)](https://doi.org/10.1021/ie502203w) -#![warn(clippy::all)] -#![allow(clippy::too_many_arguments)] #[cfg(feature = "dft")] mod dft; diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs index 437226f69..960dfc8c2 100644 --- a/src/gc_pcsaft/python/mod.rs +++ b/src/gc_pcsaft/python/mod.rs @@ -2,7 +2,6 @@ use super::dft::GcPcSaftFunctionalParameters; use super::eos::GcPcSaftEosParameters; use super::record::GcPcSaftRecord; -use crate::association::PyAssociationRecord; use feos_core::parameter::{ BinaryRecord, IdentifierOption, ParameterError, ParameterHetero, SegmentRecord, }; @@ -26,15 +25,18 @@ pub struct PyGcPcSaftRecord(GcPcSaftRecord); impl PyGcPcSaftRecord { #[new] #[pyo3( - text_signature = "(m, sigma, epsilon_k, mu=None, association_record=None, psi_dft=None)", - signature = (m, sigma, epsilon_k, mu=None, association_record=None, psi_dft=None) + text_signature = "(m, sigma, epsilon_k, mu=None, kappa_ab=None, epsilon_k_ab=None, na=None, nb=None, nc=None, psi_dft=None)" )] fn new( m: f64, sigma: f64, epsilon_k: f64, mu: Option, - association_record: Option, + kappa_ab: Option, + epsilon_k_ab: Option, + na: Option, + nb: Option, + nc: Option, psi_dft: Option, ) -> Self { Self(GcPcSaftRecord::new( @@ -42,7 +44,11 @@ impl PyGcPcSaftRecord { sigma, epsilon_k, mu, - association_record.map(|r| r.0), + kappa_ab, + epsilon_k_ab, + na, + nb, + nc, psi_dft, )) } @@ -68,8 +74,28 @@ impl PyGcPcSaftRecord { } #[getter] - fn get_association_record(&self) -> Option { - self.0.association_record.map(PyAssociationRecord) + fn get_kappa_ab(&self) -> Option { + self.0.association_record.map(|a| a.parameters.kappa_ab) + } + + #[getter] + fn get_epsilon_k_ab(&self) -> Option { + self.0.association_record.map(|a| a.parameters.epsilon_k_ab) + } + + #[getter] + fn get_na(&self) -> Option { + self.0.association_record.map(|a| a.na) + } + + #[getter] + fn get_nb(&self) -> Option { + self.0.association_record.map(|a| a.nb) + } + + #[getter] + fn get_nc(&self) -> Option { + self.0.association_record.map(|a| a.nc) } fn __repr__(&self) -> PyResult { @@ -148,7 +174,6 @@ pub fn gc_pcsaft(_py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; - m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/src/gc_pcsaft/record.rs b/src/gc_pcsaft/record.rs index 93c4fadb9..43b043d49 100644 --- a/src/gc_pcsaft/record.rs +++ b/src/gc_pcsaft/record.rs @@ -1,4 +1,5 @@ use crate::association::AssociationRecord; +use num_traits::Zero; use serde::{Deserialize, Serialize}; /// gc-PC-SAFT pure-component parameters. @@ -16,7 +17,7 @@ pub struct GcPcSaftRecord { /// Association parameters #[serde(flatten)] #[serde(skip_serializing_if = "Option::is_none")] - pub association_record: Option, + pub association_record: Option>, /// interaction range parameter for the dispersion functional #[serde(skip_serializing_if = "Option::is_none")] pub psi_dft: Option, @@ -28,9 +29,31 @@ impl GcPcSaftRecord { sigma: f64, epsilon_k: f64, mu: Option, - association_record: Option, + kappa_ab: Option, + epsilon_k_ab: Option, + na: Option, + nb: Option, + nc: Option, psi_dft: Option, ) -> Self { + let association_record = if kappa_ab.is_none() + && epsilon_k_ab.is_none() + && na.is_none() + && nb.is_none() + && nc.is_none() + { + None + } else { + Some(AssociationRecord::new( + GcPcSaftAssociationRecord::new( + kappa_ab.unwrap_or_default(), + epsilon_k_ab.unwrap_or_default(), + ), + na.unwrap_or_default(), + nb.unwrap_or_default(), + nc.unwrap_or_default(), + )) + }; Self { m, sigma, @@ -56,3 +79,31 @@ impl std::fmt::Display for GcPcSaftRecord { write!(f, ")") } } + +#[derive(Serialize, Deserialize, Clone, Copy, Default)] +pub struct GcPcSaftAssociationRecord { + /// Association volume parameter + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub kappa_ab: f64, + /// Association energy parameter in units of Kelvin + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub epsilon_k_ab: f64, +} + +impl GcPcSaftAssociationRecord { + pub fn new(kappa_ab: f64, epsilon_k_ab: f64) -> Self { + Self { + kappa_ab, + epsilon_k_ab, + } + } +} + +impl std::fmt::Display for GcPcSaftAssociationRecord { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "GcPcSaftAssociationRecord(kappa_ab={}", self.kappa_ab)?; + write!(f, ", epsilon_k_ab={})", self.epsilon_k_ab) + } +} diff --git a/src/lib.rs b/src/lib.rs index b7aeecb57..2f5659b2d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -50,12 +50,12 @@ pub mod association; pub mod hard_sphere; // models +#[cfg(feature = "epcsaft")] +pub mod epcsaft; #[cfg(feature = "gc_pcsaft")] pub mod gc_pcsaft; #[cfg(feature = "pcsaft")] pub mod pcsaft; -#[cfg(feature = "epcsaft")] -pub mod epcsaft; #[cfg(feature = "pets")] pub mod pets; #[cfg(feature = "saftvrmie")] diff --git a/src/pcsaft/mod.rs b/src/pcsaft/mod.rs index 5d3744b4f..918301d7e 100644 --- a/src/pcsaft/mod.rs +++ b/src/pcsaft/mod.rs @@ -1,8 +1,6 @@ //! Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) //! //! [Gross et al. (2001)](https://doi.org/10.1021/ie0003887) -#![warn(clippy::all)] -#![allow(clippy::too_many_arguments)] #[cfg(feature = "dft")] mod dft; diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index da72435b8..1d8899d6a 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -1,4 +1,6 @@ -use crate::association::{AssociationParameters, AssociationRecord, BinaryAssociationRecord}; +use crate::association::{ + AssociationParameters, AssociationRecord, AssociationStrength, BinaryAssociationRecord, +}; use crate::hard_sphere::{HardSphereProperties, MonomerShape}; use conv::ValueInto; use feos_core::parameter::{ @@ -11,6 +13,7 @@ use num_traits::Zero; use serde::{Deserialize, Serialize}; use std::collections::HashMap; use std::fmt::Write; +use std::sync::Arc; /// PC-SAFT pure-component parameters. #[derive(Serialize, Deserialize, Clone, Default)] @@ -30,7 +33,7 @@ pub struct PcSaftRecord { /// Association parameters #[serde(flatten)] #[serde(skip_serializing_if = "Option::is_none")] - pub association_record: Option, + pub association_record: Option>, /// Entropy scaling coefficients for the viscosity #[serde(skip_serializing_if = "Option::is_none")] pub viscosity: Option<[f64; 4]>, @@ -67,8 +70,8 @@ impl FromSegments for PcSaftRecord { .filter_map(|(s, n)| { s.association_record.as_ref().map(|record| { [ - record.kappa_ab * n, - record.epsilon_k_ab * n, + record.parameters.kappa_ab * n, + record.parameters.epsilon_k_ab * n, record.na * n, record.nb * n, record.nc * n, @@ -85,7 +88,12 @@ impl FromSegments for PcSaftRecord { ] }) .map(|[kappa_ab, epsilon_k_ab, na, nb, nc]| { - AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc) + AssociationRecord::new( + PcSaftAssociationRecord::new(kappa_ab, epsilon_k_ab), + na, + nb, + nc, + ) }); // entropy scaling @@ -239,23 +247,18 @@ impl PcSaftRecord { diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, ) -> PcSaftRecord { - let association_record = if kappa_ab.is_none() - && epsilon_k_ab.is_none() - && na.is_none() - && nb.is_none() - && nc.is_none() - { - None - } else { - Some(AssociationRecord::new( - kappa_ab.unwrap_or_default(), - epsilon_k_ab.unwrap_or_default(), - na.unwrap_or_default(), - nb.unwrap_or_default(), - nc.unwrap_or_default(), - )) - }; - PcSaftRecord { + let association_record = + if let (Some(kappa_ab), Some(epsilon_k_ab)) = (kappa_ab, epsilon_k_ab) { + Some(AssociationRecord::new( + PcSaftAssociationRecord::new(kappa_ab, epsilon_k_ab), + na.unwrap_or_default(), + nb.unwrap_or_default(), + nc.unwrap_or_default(), + )) + } else { + None + }; + Self { m, sigma, epsilon_k, @@ -269,8 +272,32 @@ impl PcSaftRecord { } } +#[derive(Serialize, Deserialize, Clone, Copy, Default)] +pub struct PcSaftAssociationRecord { + /// Association volume parameter + pub kappa_ab: f64, + /// Association energy parameter in units of Kelvin + pub epsilon_k_ab: f64, +} + +impl PcSaftAssociationRecord { + pub fn new(kappa_ab: f64, epsilon_k_ab: f64) -> Self { + Self { + kappa_ab, + epsilon_k_ab, + } + } +} + +impl std::fmt::Display for PcSaftAssociationRecord { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "PcSaftAssociationRecord(kappa_ab={}", self.kappa_ab)?; + write!(f, ", epsilon_k_ab={})", self.epsilon_k_ab) + } +} + /// PC-SAFT binary interaction parameters. -#[derive(Serialize, Deserialize, Clone, Default)] +#[derive(Serialize, Deserialize, Clone, Copy, Default)] pub struct PcSaftBinaryRecord { /// Binary dispersion interaction parameter #[serde(skip_serializing_if = "f64::is_zero")] @@ -278,7 +305,7 @@ pub struct PcSaftBinaryRecord { pub k_ij: f64, /// Binary association parameters #[serde(flatten)] - association: Option, + association: Option>, } impl From for PcSaftBinaryRecord { @@ -302,7 +329,10 @@ impl PcSaftBinaryRecord { let association = if kappa_ab.is_none() && epsilon_k_ab.is_none() { None } else { - Some(BinaryAssociationRecord::new(kappa_ab, epsilon_k_ab, None)) + Some(BinaryAssociationRecord::new( + PcSaftBinaryAssociationRecord::new(kappa_ab, epsilon_k_ab), + None, + )) }; Self { k_ij, association } } @@ -328,10 +358,10 @@ impl std::fmt::Display for PcSaftBinaryRecord { tokens.push(format!("k_ij={}", self.k_ij)); } if let Some(association) = self.association { - if let Some(kappa_ab) = association.kappa_ab { + if let Some(kappa_ab) = association.parameters.kappa_ab { tokens.push(format!("kappa_ab={}", kappa_ab)); } - if let Some(epsilon_k_ab) = association.epsilon_k_ab { + if let Some(epsilon_k_ab) = association.parameters.epsilon_k_ab { tokens.push(format!("epsilon_k_ab={}", epsilon_k_ab)); } } @@ -339,6 +369,25 @@ impl std::fmt::Display for PcSaftBinaryRecord { } } +#[derive(Serialize, Deserialize, Clone, Copy, Default)] +pub struct PcSaftBinaryAssociationRecord { + /// Cross-association association volume parameter. + #[serde(skip_serializing_if = "Option::is_none")] + pub kappa_ab: Option, + /// Cross-association energy parameter. + #[serde(skip_serializing_if = "Option::is_none")] + pub epsilon_k_ab: Option, +} + +impl PcSaftBinaryAssociationRecord { + pub fn new(kappa_ab: Option, epsilon_k_ab: Option) -> Self { + Self { + kappa_ab, + epsilon_k_ab, + } + } +} + /// Parameter set required for the PC-SAFT equation of state and Helmholtz energy functional. pub struct PcSaftParameters { pub molarweight: Array1, @@ -349,7 +398,7 @@ pub struct PcSaftParameters { pub q: Array1, pub mu2: Array1, pub q2: Array1, - pub association: AssociationParameters, + pub association: Arc>, pub sigma_ij: Array2, pub epsilon_k_ij: Array2, pub e_k_ij: Array2, @@ -425,11 +474,11 @@ impl Parameter for PcSaftParameters { .iter() .flat_map(|r| { r.indexed_iter() - .filter_map(|(i, record)| record.association.map(|r| (i, r))) + .filter_map(|((i, j), record)| record.association.map(|r| ([i, j], r))) }) .collect(); let association = - AssociationParameters::new(&association_records, &sigma, &binary_association, None); + AssociationParameters::new(&association_records, &binary_association, None); let k_ij = binary_records.as_ref().map(|br| br.map(|br| br.k_ij)); let mut sigma_ij = Array::zeros((n, n)); @@ -485,7 +534,7 @@ impl Parameter for PcSaftParameters { q, mu2, q2, - association, + association: Arc::new(association), sigma_ij, epsilon_k_ij, e_k_ij, @@ -524,6 +573,41 @@ impl HardSphereProperties for PcSaftParameters { } } +impl AssociationStrength for PcSaftParameters { + type Record = PcSaftAssociationRecord; + type BinaryRecord = PcSaftBinaryAssociationRecord; + + fn association_strength + Copy>( + &self, + temperature: D, + comp_i: usize, + comp_j: usize, + assoc_ij: Self::Record, + ) -> D { + let si = self.sigma[comp_i]; + let sj = self.sigma[comp_j]; + (temperature.recip() * assoc_ij.epsilon_k_ab).exp_m1() + * assoc_ij.kappa_ab + * (si * sj).powf(1.5) + } + + fn combining_rule(parameters_i: Self::Record, parameters_j: Self::Record) -> Self::Record { + Self::Record { + kappa_ab: (parameters_i.kappa_ab * parameters_j.kappa_ab).sqrt(), + epsilon_k_ab: 0.5 * (parameters_i.epsilon_k_ab + parameters_j.epsilon_k_ab), + } + } + + fn update_binary(parameters_ij: &mut Self::Record, binary_parameters: Self::BinaryRecord) { + if let Some(kappa_ab) = binary_parameters.kappa_ab { + parameters_ij.kappa_ab = kappa_ab + } + if let Some(epsilon_k_ab) = binary_parameters.epsilon_k_ab { + parameters_ij.epsilon_k_ab = epsilon_k_ab + } + } +} + impl PcSaftParameters { pub fn to_markdown(&self) -> String { let mut output = String::new(); @@ -536,10 +620,9 @@ impl PcSaftParameters { for (i, record) in self.pure_records.iter().enumerate() { let component = record.identifier.name.clone(); let component = component.unwrap_or(format!("Component {}", i + 1)); - let association = record - .model_record - .association_record - .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, 0.0, 0.0, 0.0)); + let association = record.model_record.association_record.unwrap_or_else(|| { + AssociationRecord::new(PcSaftAssociationRecord::new(0.0, 0.0), 0.0, 0.0, 0.0) + }); write!( o, "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", @@ -550,8 +633,8 @@ impl PcSaftParameters { record.model_record.epsilon_k, record.model_record.mu.unwrap_or(0.0), record.model_record.q.unwrap_or(0.0), - association.kappa_ab, - association.epsilon_k_ab, + association.parameters.kappa_ab, + association.parameters.epsilon_k_ab, association.na, association.nb, association.nc diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index a65efaa71..211c4823d 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -111,12 +111,12 @@ impl PyPcSaftRecord { #[getter] fn get_kappa_ab(&self) -> Option { - self.0.association_record.map(|a| a.kappa_ab) + self.0.association_record.map(|a| a.parameters.kappa_ab) } #[getter] fn get_epsilon_k_ab(&self) -> Option { - self.0.association_record.map(|a| a.epsilon_k_ab) + self.0.association_record.map(|a| a.parameters.epsilon_k_ab) } #[getter] diff --git a/src/pets/mod.rs b/src/pets/mod.rs index 5f496b9b4..92e0f1242 100644 --- a/src/pets/mod.rs +++ b/src/pets/mod.rs @@ -5,10 +5,6 @@ //! PeTS is an equation of state for the truncated and shifted Lennar-Jones fluid with cut-off //! distance 2.5 $\sigma$. //! It utilizes a hard-sphere fluid as reference with an attactive perturbation. -#![warn(clippy::all)] -#![allow(clippy::too_many_arguments)] -#![allow(clippy::many_single_char_names)] -#![allow(clippy::suspicious_operation_groupings)] #[cfg(feature = "dft")] mod dft; diff --git a/src/python/eos.rs b/src/python/eos.rs index e7f3acae5..422323deb 100644 --- a/src/python/eos.rs +++ b/src/python/eos.rs @@ -1,4 +1,8 @@ use crate::eos::ResidualModel; +#[cfg(feature = "epcsaft")] +use crate::epcsaft::python::PyElectrolytePcSaftParameters; +#[cfg(feature = "epcsaft")] +use crate::epcsaft::{ElectrolytePcSaft, ElectrolytePcSaftOptions, ElectrolytePcSaftVariants}; #[cfg(feature = "estimator")] use crate::estimator::*; #[cfg(feature = "gc_pcsaft")] @@ -14,10 +18,6 @@ use crate::impl_estimator_entropy_scaling; use crate::pcsaft::python::PyPcSaftParameters; #[cfg(feature = "pcsaft")] use crate::pcsaft::{DQVariants, PcSaft, PcSaftOptions}; -#[cfg(feature = "epcsaft")] -use crate::epcsaft::python::PyElectrolytePcSaftParameters; -#[cfg(feature = "epcsaft")] -use crate::epcsaft::{ElectrolytePcSaft, ElectrolytePcSaftOptions, ElectrolytePcSaftVariants}; #[cfg(feature = "pets")] use crate::pets::python::PyPetsParameters; #[cfg(feature = "pets")] @@ -212,20 +212,20 @@ impl PyEquationOfState { /// Returns /// ------- /// EquationOfState - /// The PC-SAFT equation of state that can be used to compute thermodynamic + /// The ePC-SAFT equation of state that can be used to compute thermodynamic /// states. #[cfg(feature = "epcsaft")] #[staticmethod] #[pyo3( signature = (parameters, max_eta=0.5, max_iter_cross_assoc=50, tol_cross_assoc=1e-10, epcsaft_variant=ElectrolytePcSaftVariants::Advanced), - text_signature = "(parameters, max_eta=0.5, max_iter_cross_assoc=50, tol_cross_assoc=1e-10, epcsaft_variant=advanced)", + text_signature = "(parameters, max_eta=0.5, max_iter_cross_assoc=50, tol_cross_assoc=1e-10, epcsaft_variant)", )] pub fn epcsaft( parameters: PyElectrolytePcSaftParameters, max_eta: f64, max_iter_cross_assoc: usize, tol_cross_assoc: f64, - epcsaft_variant: ElectrolytePcSaftVariants + epcsaft_variant: ElectrolytePcSaftVariants, ) -> Self { let options = ElectrolytePcSaftOptions { max_eta, @@ -233,10 +233,9 @@ impl PyEquationOfState { tol_cross_assoc, epcsaft_variant, }; - let residual = Arc::new(ResidualModel::ElectrolytePcSaft(ElectrolytePcSaft::with_options( - parameters.0, - options, - ))); + let residual = Arc::new(ResidualModel::ElectrolytePcSaft( + ElectrolytePcSaft::with_options(parameters.0, options), + )); let ideal_gas = Arc::new(IdealGasModel::NoModel(residual.components())); Self(Arc::new(EquationOfState::new(ideal_gas, residual))) } diff --git a/src/saftvrmie/python.rs b/src/saftvrmie/python.rs index e5a830cb2..a48c6d682 100644 --- a/src/saftvrmie/python.rs +++ b/src/saftvrmie/python.rs @@ -262,8 +262,6 @@ impl PySaftVRMieParameters { pub fn saftvrmie(_py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; - m.add_class::()?; - m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/src/saftvrqmie/mod.rs b/src/saftvrqmie/mod.rs index 158094933..c6f89e521 100644 --- a/src/saftvrqmie/mod.rs +++ b/src/saftvrqmie/mod.rs @@ -7,8 +7,7 @@ //! # Literature //! - Pure substances: [Aasen et al. (2019)](https://aip.scitation.org/doi/10.1063/1.5111364) //! - Binary mixtures: [Aasen et al. (2020)](https://aip.scitation.org/doi/10.1063/1.5136079) -#![warn(clippy::all)] -#![allow(clippy::too_many_arguments)] + #[cfg(feature = "dft")] mod dft; mod eos;