Allow multiple association sites per molecule in PC-SAFT

feos-org · prehner · Jun 20, 2025 · Apr 25, 2025 · Apr 30, 2025 · Apr 30, 2025
commit 4bbe0bf53742842d7ac4713da2652aca38a48194
diff --git a/Cargo.toml b/Cargo.toml
@@ -42,6 +42,7 @@ libm = "0.2"
 gauss-quad = "0.2"
 approx = "0.5"
 criterion = "0.5"
+arrayvec = "0.7"
 
 feos-core = { version = "0.8", path = "crates/feos-core" }
 feos-dft = { version = "0.8", path = "crates/feos-dft" }

diff --git a/crates/feos-core/Cargo.toml b/crates/feos-core/Cargo.toml
@@ -26,6 +26,7 @@ serde_json = { workspace = true, features = ["preserve_order"] }
 indexmap = { workspace = true, features = ["serde"] }
 rayon = { workspace = true, optional = true }
 typenum = { workspace = true }
+itertools = { workspace = true }
 
 [dev-dependencies]
 approx = { workspace = true }

diff --git a/crates/feos-core/src/cubic.rs b/crates/feos-core/src/cubic.rs
@@ -62,6 +62,8 @@ pub struct PengRobinsonParameters {
     molarweight: Array1<f64>,
     /// List of pure component records
     pure_records: Vec<PureRecord<PengRobinsonRecord>>,
+    /// List of binary records
+    binary_records: Vec<([usize; 2], f64)>,
 }
 
 impl std::fmt::Display for PengRobinsonParameters {
@@ -100,7 +102,7 @@ impl PengRobinsonParameters {
                 PureRecord::new(id, molarweight[i], record)
             })
             .collect();
-        PengRobinsonParameters::from_records(records, None)
+        PengRobinsonParameters::from_records(records, vec![])
     }
 }
 
@@ -111,7 +113,7 @@ impl Parameter for PengRobinsonParameters {
     /// Creates parameters from pure component records.
     fn from_records(
         pure_records: Vec<PureRecord<Self::Pure>>,
-        binary_records: Option<Array2<Self::Binary>>,
+        binary_records: Vec<([usize; 2], Self::Binary)>,
     ) -> FeosResult<Self> {
         let n = pure_records.len();
 
@@ -130,7 +132,11 @@ impl Parameter for PengRobinsonParameters {
             kappa[i] = 0.37464 + (1.54226 - 0.26992 * r.acentric_factor) * r.acentric_factor;
         }
 
-        let k_ij = binary_records.unwrap_or_else(|| Array2::zeros([n; 2]));
+        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;
+        }
 
         Ok(Self {
             tc,
@@ -140,11 +146,12 @@ impl Parameter for PengRobinsonParameters {
             kappa,
             molarweight,
             pure_records,
+            binary_records,
         })
     }
 
-    fn records(&self) -> (&[PureRecord<PengRobinsonRecord>], Option<&Array2<f64>>) {
-        (&self.pure_records, Some(&self.k_ij))
+    fn records(&self) -> (&[PureRecord<PengRobinsonRecord>], &[([usize; 2], f64)]) {
+        (&self.pure_records, &self.binary_records)
     }
 }
 

diff --git a/crates/feos-core/src/parameter/mod.rs b/crates/feos-core/src/parameter/mod.rs
@@ -2,7 +2,7 @@
 
 use crate::errors::*;
 use indexmap::{IndexMap, IndexSet};
-use ndarray::Array2;
+use itertools::Itertools;
 use serde::de::DeserializeOwned;
 use serde::{Deserialize, Serialize};
 use std::collections::HashMap;
@@ -36,12 +36,12 @@ where
     /// Creates parameters from records for pure substances and possibly binary parameters.
     fn from_records(
         pure_records: Vec<PureRecord<Self::Pure>>,
-        binary_records: Option<Array2<Self::Binary>>,
+        binary_records: Vec<([usize; 2], Self::Binary)>,
     ) -> FeosResult<Self>;
 
     /// Creates parameters for a pure component from a pure record.
     fn new_pure(pure_record: PureRecord<Self::Pure>) -> FeosResult<Self> {
-        Self::from_records(vec![pure_record], None)
+        Self::from_records(vec![pure_record], vec![])
     }
 
     /// Creates parameters for a binary system from pure records and an optional
@@ -50,15 +50,7 @@ where
         pure_records: Vec<PureRecord<Self::Pure>>,
         binary_record: Option<Self::Binary>,
     ) -> FeosResult<Self> {
-        let binary_record = binary_record.map(|br| {
-            Array2::from_shape_fn([2, 2], |(i, j)| {
-                if i == j {
-                    Self::Binary::default()
-                } else {
-                    br.clone()
-                }
-            })
-        });
+        let binary_record = binary_record.into_iter().map(|r| ([0, 0], r)).collect();
         Self::from_records(pure_records, binary_record)
     }
 
@@ -69,27 +61,22 @@ where
             .into_iter()
             .map(|r| PureRecord::new(Default::default(), Default::default(), r))
             .collect();
-        Self::from_records(pure_records, None)
+        Self::from_records(pure_records, vec![])
     }
 
     /// Return the original pure and binary records that were used to construct the parameters.
     #[expect(clippy::type_complexity)]
-    fn records(&self) -> (&[PureRecord<Self::Pure>], Option<&Array2<Self::Binary>>);
+    fn records(&self) -> (&[PureRecord<Self::Pure>], &[([usize; 2], Self::Binary)]);
 
     /// Helper function to build matrix from list of records in correct order.
     ///
     /// If the identifiers in `binary_records` are not a subset of those in
     /// `pure_records`, the `Default` implementation of Self::Binary is used.
-    #[expect(clippy::expect_fun_call)]
     fn binary_matrix_from_records(
         pure_records: &[PureRecord<Self::Pure>],
         binary_records: &[BinaryRecord<Self::Binary>],
         identifier_option: IdentifierOption,
-    ) -> Option<Array2<Self::Binary>> {
-        if binary_records.is_empty() {
-            return None;
-        }
-
+    ) -> Vec<([usize; 2], Self::Binary)> {
         // Build Hashmap (id, id) -> BinaryRecord
         let binary_map: HashMap<_, _> = {
             binary_records
@@ -101,28 +88,31 @@ where
                 })
                 .collect()
         };
-        let n = pure_records.len();
-        Some(Array2::from_shape_fn([n, n], |(i, j)| {
-            let id1 = pure_records[i]
-                .identifier
-                .as_str(identifier_option)
-                .expect(&format!(
-                    "No identifier for given identifier_option for pure record {}.",
-                    i
-                ));
-            let id2 = pure_records[j]
-                .identifier
-                .as_str(identifier_option)
-                .expect(&format!(
-                    "No identifier for given identifier_option for pure record {}.",
-                    j
-                ));
-            binary_map
-                .get(&(id1, id2))
-                .or_else(|| binary_map.get(&(id2, id1)))
-                .cloned()
-                .unwrap_or_default()
-        }))
+
+        // look up pure records in Hashmap
+        pure_records
+            .iter()
+            .enumerate()
+            .array_combinations()
+            .filter_map(|[(i1, p1), (i2, p2)]| {
+                let id1 = p1.identifier.as_str(identifier_option).unwrap_or_else(|| {
+                    panic!(
+                        "No {} for pure record {} ({}).",
+                        identifier_option, i1, p1.identifier
+                    )
+                });
+                let id2 = p2.identifier.as_str(identifier_option).unwrap_or_else(|| {
+                    panic!(
+                        "No {} for pure record {} ({}).",
+                        identifier_option, i2, p2.identifier
+                    )
+                });
+                binary_map
+                    .get(&(id1, id2))
+                    .or_else(|| binary_map.get(&(id2, id1)))
+                    .map(|br| ([i1, i2], br.clone()))
+            })
+            .collect()
     }
 
     /// Creates parameters from substance information stored in json files.
@@ -230,13 +220,14 @@ where
         // full matrix of binary records from the gc method.
         // If a specific segment-segment interaction is not in the binary map,
         // the default value is used.
-        let n = pure_records.len();
-        let mut binary_records = Array2::default([n, n]);
-        for i in 0..n {
-            for j in i + 1..n {
+        let x = segment_counts
+            .iter()
+            .enumerate()
+            .array_combinations()
+            .map(|[(i, sc1), (j, sc2)]| {
                 let mut vec = Vec::new();
-                for (id1, &n1) in segment_counts[i].iter() {
-                    for (id2, &n2) in segment_counts[j].iter() {
+                for (id1, &n1) in sc1.iter() {
+                    for (id2, &n2) in sc2.iter() {
                         let binary = binary_map
                             .get(&(id1.clone(), id2.clone()))
                             .or_else(|| binary_map.get(&(id2.clone(), id1.clone())))
@@ -245,13 +236,11 @@ where
                         vec.push((binary, n1, n2));
                     }
                 }
-                let kij = Self::Binary::from_segments_binary(&vec)?;
-                binary_records[(i, j)] = kij.clone();
-                binary_records[(j, i)] = kij;
-            }
-        }
+                Self::Binary::from_segments_binary(&vec).map(|kij| ([i, j], kij))
+            })
+            .collect::<Result<_, _>>()?;
 
-        Self::from_records(pure_records, Some(binary_records))
+        Self::from_records(pure_records, x)
     }
 
     /// Creates parameters from segment information stored in json files.
@@ -332,12 +321,9 @@ where
             .iter()
             .map(|&i| pure_records[i].clone())
             .collect();
-        let n = component_list.len();
-        let binary_records = binary_records.map(|br| {
-            Array2::from_shape_fn([n, n], |(i, j)| {
-                br[(component_list[i], component_list[j])].clone()
-            })
-        });
+        let mut binary_records = binary_records.to_vec();
+        binary_records
+            .retain(|([i, j], _)| component_list.contains(i) && component_list.contains(j));
 
         Self::from_records(pure_records, binary_records)
             .expect("failed to create subset from parameters.")

diff --git a/crates/feos-core/src/parameter/model_record.rs b/crates/feos-core/src/parameter/model_record.rs
@@ -108,7 +108,7 @@ where
         write!(f, "PureRecord(")?;
         write!(f, "\n\tidentifier={},", self.identifier)?;
         write!(f, "\n\tmolarweight={},", self.molarweight)?;
-        write!(f, "\n\tmodel_record={},", self.model_record)?;
+        write!(f, "\n\tmodel_record={}", self.model_record)?;
         write!(f, "\n)")
     }
 }

diff --git a/crates/feos-core/tests/parameters.rs b/crates/feos-core/tests/parameters.rs
@@ -1,39 +1,38 @@
 use feos_core::FeosError;
 use feos_core::FeosResult;
 use feos_core::parameter::*;
-use ndarray::Array2;
 use serde::{Deserialize, Serialize};
 
 #[derive(Debug, Clone, Serialize, Deserialize, Default)]
 struct MyPureModel {
     a: f64,
 }
 
-#[derive(Debug, Clone, Serialize, Deserialize, Default, PartialEq)]
+#[derive(Debug, Clone, Copy, Serialize, Deserialize, Default, PartialEq)]
 struct MyBinaryModel {
     b: f64,
 }
 
 struct MyParameter {
     pure_records: Vec<PureRecord<MyPureModel>>,
-    binary_records: Option<Array2<MyBinaryModel>>,
+    binary_records: Vec<([usize; 2], MyBinaryModel)>,
 }
 
 impl Parameter for MyParameter {
     type Pure = MyPureModel;
     type Binary = MyBinaryModel;
     fn from_records(
         pure_records: Vec<PureRecord<MyPureModel>>,
-        binary_records: Option<Array2<MyBinaryModel>>,
+        binary_records: Vec<([usize; 2], MyBinaryModel)>,
     ) -> FeosResult<Self> {
         Ok(Self {
             pure_records,
             binary_records,
         })
     }
 
-    fn records(&self) -> (&[PureRecord<MyPureModel>], Option<&Array2<MyBinaryModel>>) {
-        (&self.pure_records, self.binary_records.as_ref())
+    fn records(&self) -> (&[PureRecord<MyPureModel>], &[([usize; 2], MyBinaryModel)]) {
+        (&self.pure_records, &self.binary_records)
     }
 }
 
@@ -87,7 +86,7 @@ fn from_records() -> FeosResult<()> {
 
     assert_eq!(p.pure_records[0].identifier.cas, Some("123-4-5".into()));
     assert_eq!(p.pure_records[1].identifier.cas, Some("678-9-1".into()));
-    assert_eq!(p.binary_records.unwrap()[[0, 1]].b, 12.0);
+    assert_eq!(p.binary_records[0].1.b, 12.0);
     Ok(())
 }
 
@@ -151,9 +150,8 @@ fn from_multiple_json_files() {
     // test_parameters1: a = 0.5
     // test_parameters2: a = 0.1 or 0.3
     assert_eq!(p.pure_records[1].model_record.a, 0.5);
-    let br = p.binary_records.as_ref().unwrap();
-    assert_eq!(br[[0, 1]].b, 12.0);
-    assert_eq!(br[[1, 0]].b, 12.0);
+    let br = p.binary_records;
+    assert_eq!(br[0].1.b, 12.0);
 }
 
 #[test]
@@ -206,9 +204,8 @@ fn from_records_missing_binary() -> FeosResult<()> {
 
     assert_eq!(p.pure_records[0].identifier.cas, Some("123-4-5".into()));
     assert_eq!(p.pure_records[1].identifier.cas, Some("678-9-1".into()));
-    let br = p.binary_records.as_ref().unwrap();
-    assert_eq!(br[[0, 1]], MyBinaryModel::default());
-    assert_eq!(br[[0, 1]].b, 0.0);
+    let br = p.binary_records;
+    assert_eq!(br.len(), 0);
     Ok(())
 }
 
@@ -272,12 +269,9 @@ fn from_records_correct_binary_order() -> FeosResult<()> {
     assert_eq!(p.pure_records[0].identifier.cas, Some("000-0-0".into()));
     assert_eq!(p.pure_records[1].identifier.cas, Some("123-4-5".into()));
     assert_eq!(p.pure_records[2].identifier.cas, Some("678-9-1".into()));
-    let br = p.binary_records.as_ref().unwrap();
-    assert_eq!(br[[0, 1]], MyBinaryModel::default());
-    assert_eq!(br[[1, 0]], MyBinaryModel::default());
-    assert_eq!(br[[0, 2]], MyBinaryModel::default());
-    assert_eq!(br[[2, 0]], MyBinaryModel::default());
-    assert_eq!(br[[2, 1]].b, 12.0);
-    assert_eq!(br[[1, 2]].b, 12.0);
+    let ([i, j], br) = p.binary_records[0];
+    assert_eq!(i, 1);
+    assert_eq!(j, 2);
+    assert_eq!(br.b, 12.0);
     Ok(())
 }
diff --git a/crates/feos/Cargo.toml b/crates/feos/Cargo.toml
@@ -22,6 +22,7 @@ indexmap = { workspace = true }
 rayon = { workspace = true, optional = true }
 itertools = { workspace = true }
 typenum = { workspace = true }
+arrayvec = { workspace = true, features = ["serde"] }
 
 feos-core = { workspace = true }
 feos-derive = { workspace = true }

diff --git a/crates/feos/benches/dual_numbers.rs b/crates/feos/benches/dual_numbers.rs
@@ -4,11 +4,9 @@
 //! on the dual number type used without the overhead of the `State`
 //! creation.
 use criterion::{Criterion, criterion_group, criterion_main};
-use feos::core::{
-    Derivative, Residual, State, StateHD,
-    parameter::{IdentifierOption, Parameter},
-};
-use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos::core::parameter::{IdentifierOption, Parameter};
+use feos::core::{Derivative, Residual, State, StateHD};
+use feos::pcsaft::{PcSaft, PcSaftBinaryRecord, PcSaftParameters};
 use ndarray::{Array, ScalarOperand, arr1};
 use num_dual::DualNum;
 use quantity::*;
@@ -115,8 +113,8 @@ fn methane_co2_pcsaft(c: &mut Criterion) {
     )
     .unwrap();
     let k_ij = -0.0192211646;
-    let parameters =
-        PcSaftParameters::new_binary(parameters.pure_records, Some(k_ij.into())).unwrap();
+    let br = PcSaftBinaryRecord::new(k_ij, None, None, vec![]);
+    let parameters = PcSaftParameters::new_binary(parameters.pure_records, Some(br)).unwrap();
     let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
 
     // 230 K, 50 bar, x0 = 0.15