Handle association outside of individual models

feos-org · prehner · Jun 20, 2025 · Apr 25, 2025 · Apr 30, 2025 · Apr 30, 2025
commit 692a754ccb2753c37642aed7b33c15206aab64f2
diff --git a/crates/feos-core/Cargo.toml b/crates/feos-core/Cargo.toml
@@ -27,6 +27,7 @@ indexmap = { workspace = true, features = ["serde"] }
 rayon = { workspace = true, optional = true }
 typenum = { workspace = true }
 itertools = { workspace = true }
+arrayvec = { workspace = true, features = ["serde"] }
 
 [dev-dependencies]
 approx = { workspace = true }

diff --git a/crates/feos-core/src/cubic.rs b/crates/feos-core/src/cubic.rs
@@ -7,7 +7,7 @@
 use crate::FeosError;
 use crate::equation_of_state::{Components, Molarweight, Residual};
 use crate::errors::FeosResult;
-use crate::parameter::{Identifier, Parameter, PureRecord};
+use crate::parameter::{BinaryRecord, Identifier, Parameter, PureRecord};
 use crate::state::StateHD;
 use ndarray::{Array1, Array2, ScalarOperand};
 use num_dual::DualNum;
@@ -61,9 +61,9 @@ pub struct PengRobinsonParameters {
     /// Molar weight in units of g/mol
     molarweight: Array1<f64>,
     /// List of pure component records
-    pure_records: Vec<PureRecord<PengRobinsonRecord>>,
+    pure_records: Vec<PureRecord<PengRobinsonRecord, ()>>,
     /// List of binary records
-    binary_records: Vec<([usize; 2], f64)>,
+    binary_records: Vec<BinaryRecord<usize, f64, ()>>,
 }
 
 impl std::fmt::Display for PengRobinsonParameters {
@@ -109,11 +109,12 @@ impl PengRobinsonParameters {
 impl Parameter for PengRobinsonParameters {
     type Pure = PengRobinsonRecord;
     type Binary = f64;
+    type Association = ();
 
     /// Creates parameters from pure component records.
     fn from_records(
-        pure_records: Vec<PureRecord<Self::Pure>>,
-        binary_records: Vec<([usize; 2], Self::Binary)>,
+        pure_records: Vec<PureRecord<Self::Pure, ()>>,
+        binary_records: Vec<BinaryRecord<usize, Self::Binary, ()>>,
     ) -> FeosResult<Self> {
         let n = pure_records.len();
 
@@ -133,9 +134,9 @@ impl Parameter for PengRobinsonParameters {
         }
 
         let mut k_ij = Array2::zeros([n; 2]);
-        for &([i, j], r) in &binary_records {
-            k_ij[[i, j]] = r;
-            k_ij[[j, i]] = r;
+        for r in &binary_records {
+            k_ij[[r.id1, r.id2]] = r.model_record.unwrap_or_default();
+            k_ij[[r.id2, r.id1]] = r.model_record.unwrap_or_default();
         }
 
         Ok(Self {
@@ -150,7 +151,12 @@ impl Parameter for PengRobinsonParameters {
         })
     }
 
-    fn records(&self) -> (&[PureRecord<PengRobinsonRecord>], &[([usize; 2], f64)]) {
+    fn records(
+        &self,
+    ) -> (
+        &[PureRecord<PengRobinsonRecord, ()>],
+        &[BinaryRecord<usize, f64, ()>],
+    ) {
         (&self.pure_records, &self.binary_records)
     }
 }
@@ -240,7 +246,7 @@ mod tests {
     use quantity::{KELVIN, PASCAL};
     use std::sync::Arc;
 
-    fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord>> {
+    fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord, ()>> {
         let records = r#"[
             {
                 "identifier": {
@@ -251,11 +257,9 @@ mod tests {
                     "inchi": "InChI=1/C3H8/c1-3-2/h3H2,1-2H3",
                     "formula": "C3H8"
                 },
-                "model_record": {
-                    "tc": 369.96,
-                    "pc": 4250000.0,
-                    "acentric_factor": 0.153
-                },
+                "tc": 369.96,
+                "pc": 4250000.0,
+                "acentric_factor": 0.153,
                 "molarweight": 44.0962
             },
             {
@@ -267,11 +271,9 @@ mod tests {
                     "inchi": "InChI=1/C4H10/c1-3-4-2/h3-4H2,1-2H3",
                     "formula": "C4H10"
                 },
-                "model_record": {
-                    "tc": 425.2,
-                    "pc": 3800000.0,
-                    "acentric_factor": 0.199
-                },
+                "tc": 425.2,
+                "pc": 3800000.0,
+                "acentric_factor": 0.199,
                 "molarweight": 58.123
             }
         ]"#;

diff --git a/crates/feos-core/src/lib.rs b/crates/feos-core/src/lib.rs
@@ -43,7 +43,6 @@ pub use state::{
     Contributions, DensityInitialization, Derivative, State, StateBuilder, StateHD, StateVec,
 };
 
-
 /// Level of detail in the iteration output.
 #[derive(Copy, Clone, PartialOrd, PartialEq, Eq)]
 pub enum Verbosity {
@@ -167,16 +166,8 @@ pub trait ReferenceSystem<
 }
 
 /// Conversion to and from reduced units
-impl<
-        Inner,
-        T: Integer,
-        L: Integer,
-        M: Integer,
-        I: Integer,
-        THETA: Integer,
-        N: Integer,
-        J: Integer,
-    > ReferenceSystem<Inner, T, L, M, I, THETA, N, J>
+impl<Inner, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    ReferenceSystem<Inner, T, L, M, I, THETA, N, J>
     for Quantity<Inner, SIUnit<T, L, M, I, THETA, N, J>>
 {
     fn from_reduced(value: Inner) -> Self
@@ -203,12 +194,12 @@ impl<
 
 #[cfg(test)]
 mod tests {
-    use crate::cubic::*;
-    use crate::equation_of_state::{Components, EquationOfState, IdealGas};
-    use crate::parameter::*;
     use crate::Contributions;
     use crate::FeosResult;
     use crate::StateBuilder;
+    use crate::cubic::*;
+    use crate::equation_of_state::{Components, EquationOfState, IdealGas};
+    use crate::parameter::*;
     use approx::*;
     use ndarray::Array1;
     use num_dual::DualNum;
@@ -238,7 +229,7 @@ mod tests {
         }
     }
 
-    fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord>> {
+    fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord, ()>> {
         let records = r#"[
             {
                 "identifier": {
@@ -249,11 +240,9 @@ mod tests {
                     "inchi": "InChI=1/C3H8/c1-3-2/h3H2,1-2H3",
                     "formula": "C3H8"
                 },
-                "model_record": {
-                    "tc": 369.96,
-                    "pc": 4250000.0,
-                    "acentric_factor": 0.153
-                },
+                "tc": 369.96,
+                "pc": 4250000.0,
+                "acentric_factor": 0.153,
                 "molarweight": 44.0962
             },
             {
@@ -265,11 +254,9 @@ mod tests {
                     "inchi": "InChI=1/C4H10/c1-3-4-2/h3-4H2,1-2H3",
                     "formula": "C4H10"
                 },
-                "model_record": {
-                    "tc": 425.2,
-                    "pc": 3800000.0,
-                    "acentric_factor": 0.199
-                },
+                "tc": 425.2,
+                "pc": 3800000.0,
+                "acentric_factor": 0.199,
                 "molarweight": 58.123
             }
         ]"#;

diff --git a/crates/feos-core/src/parameter/association.rs b/crates/feos-core/src/parameter/association.rs
@@ -0,0 +1,73 @@
+use arrayvec::ArrayString;
+use num_traits::Zero;
+use serde::{Deserialize, Serialize};
+
+type SiteId = ArrayString<8>;
+
+/// Pure component association parameters.
+#[derive(Serialize, Deserialize, Clone, Copy)]
+pub struct AssociationRecord<A> {
+    #[serde(skip_serializing_if = "SiteId::is_empty")]
+    #[serde(default)]
+    pub id: SiteId,
+    #[serde(flatten)]
+    pub parameters: Option<A>,
+    /// \# of association sites of type A
+    #[serde(skip_serializing_if = "f64::is_zero")]
+    #[serde(default)]
+    pub na: f64,
+    /// \# of association sites of type B
+    #[serde(skip_serializing_if = "f64::is_zero")]
+    #[serde(default)]
+    pub nb: f64,
+    /// \# of association sites of type C
+    #[serde(skip_serializing_if = "f64::is_zero")]
+    #[serde(default)]
+    pub nc: f64,
+}
+
+impl<A> AssociationRecord<A> {
+    pub fn new(parameters: Option<A>, na: f64, nb: f64, nc: f64) -> Self {
+        Self::with_id(Default::default(), parameters, na, nb, nc)
+    }
+
+    pub fn with_id(id: SiteId, parameters: Option<A>, na: f64, nb: f64, nc: f64) -> Self {
+        Self {
+            id,
+            parameters,
+            na,
+            nb,
+            nc,
+        }
+    }
+}
+
+/// Binary association parameters.
+#[derive(Serialize, Deserialize, Clone, Copy)]
+pub struct BinaryAssociationRecord<A> {
+    // Identifier of the association site on the first molecule.
+    #[serde(skip_serializing_if = "SiteId::is_empty")]
+    #[serde(default)]
+    pub id1: SiteId,
+    // Identifier of the association site on the second molecule.
+    #[serde(skip_serializing_if = "SiteId::is_empty")]
+    #[serde(default)]
+    pub id2: SiteId,
+    // Binary association parameters
+    #[serde(flatten)]
+    pub parameters: A,
+}
+
+impl<A> BinaryAssociationRecord<A> {
+    pub fn new(parameters: A) -> Self {
+        Self::with_id(Default::default(), Default::default(), parameters)
+    }
+
+    pub fn with_id(id1: SiteId, id2: SiteId, parameters: A) -> Self {
+        Self {
+            id1,
+            id2,
+            parameters,
+        }
+    }
+}
diff --git a/crates/feos-core/src/parameter/chemical_record.rs b/crates/feos-core/src/parameter/chemical_record.rs
@@ -1,12 +1,9 @@
-use super::identifier::Identifier;
-use super::segment::SegmentRecord;
+use super::{Identifier, SegmentRecord};
 use crate::{FeosError, FeosResult};
 use num_traits::NumAssign;
 use serde::{Deserialize, Serialize};
-use std::{
-    borrow::Cow,
-    collections::{HashMap, HashSet},
-};
+use std::borrow::Cow;
+use std::collections::{HashMap, HashSet};
 
 // Auxiliary structure used to deserialize chemical records without explicit bond information.
 #[derive(Serialize, Deserialize)]
@@ -135,19 +132,20 @@ pub trait SegmentCount {
     /// molecule.
     ///
     /// The map contains the segment record as key and the count as value.
-    fn segment_map<M: Clone>(
+    #[expect(clippy::type_complexity)]
+    fn segment_map<M: Clone, A: Clone>(
         &self,
-        segment_records: &[SegmentRecord<M>],
-    ) -> FeosResult<HashMap<SegmentRecord<M>, Self::Count>> {
+        segment_records: &[SegmentRecord<M, A>],
+    ) -> FeosResult<Vec<(SegmentRecord<M, A>, Self::Count)>> {
         let count = self.segment_count();
-        let queried: HashSet<_> = count.keys().cloned().collect();
-        let mut segments: HashMap<String, SegmentRecord<M>> = segment_records
+        let queried: HashSet<_> = count.keys().collect();
+        let mut segments: HashMap<_, SegmentRecord<M, A>> = segment_records
             .iter()
-            .map(|r| (r.identifier.clone(), r.clone()))
+            .map(|r| (&r.identifier, r.clone()))
             .collect();
-        let available = segments.keys().cloned().collect();
+        let available = segments.keys().copied().collect();
         if !queried.is_subset(&available) {
-            let missing: Vec<String> = queried.difference(&available).cloned().collect();
+            let missing: Vec<_> = queried.difference(&available).collect();
             let msg = format!("{:?}", missing);
             return Err(FeosError::ComponentsNotFound(msg));
         };