Group contribution method for binary interaction parameters

feos-org · prehner · Jan 27, 2023 · Jan 7, 2023 · Jan 7, 2023 · Jan 7, 2023
commit bf56d5c328ad686f9581396ead0fd2b7c1ab070e
diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md
@@ -10,7 +10,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `StateVec::moles` getter. [#113](https://github.com/feos-org/feos/pull/113)
 - Added public constructors `PhaseDiagram::new` and `StateVec::new` that allow the creation of the respective structs from a list of `PhaseEquilibrium`s or `State`s in Rust and Python. [#113](https://github.com/feos-org/feos/pull/113)
 - Added new variant `EosError::Error` which allows dispatching generic errors that are not covered by the existing variants. [#98](https://github.com/feos-org/feos/pull/98)
-- Removed generics for units in all structs and traits in favor of static SI units. [#115](https://github.com/feos-org/feos/pull/115)
+- Added `binary_records` getter for parameter classes in Python. [#104](https://github.com/feos-org/feos/pull/104)
+- Added `BinaryRecord::from_json` and `BinarySegmentRecord::from_json` that read a list of records from a file. [#104](https://github.com/feos-org/feos/pull/104)
 
 ### Changed
 - Added `Sync` and `Send` as supertraits to `EquationOfState`. [#57](https://github.com/feos-org/feos/pull/57)
@@ -22,12 +23,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - The critical point algorithm now uses vector dual numbers to reduce the number of model evaluations and computation times. [#96](https://github.com/feos-org/feos/pull/96)
 - Renamed `State::molar_volume` to `State::partial_molar_volume` and `State::ln_phi_pure` to `State::ln_phi_pure_liquid`. [#107](https://github.com/feos-org/feos/pull/107)
 - Added a const generic parameter to `PhaseDiagram` that accounts for the number of phases analogously to `PhaseEquilibrium`. [#113](https://github.com/feos-org/feos/pull/113)
+- Removed generics for units in all structs and traits in favor of static SI units. [#115](https://github.com/feos-org/feos/pull/115)
 
 ### Packaging
 - Updated `pyo3` and `numpy` dependencies to 0.18. [#119](https://github.com/feos-org/feos/pull/119)
 - Updated `quantity` dependency to 0.6. [#119](https://github.com/feos-org/feos/pull/119)
 - Updated `num-dual` dependency to 0.6. [#119](https://github.com/feos-org/feos/pull/119)
 
+### Fixed
+- `Parameter::from_segments` only calculates binary interaction parameters for unlike molecules. [#104](https://github.com/feos-org/feos/pull/104)
+
 ## [0.3.1] - 2022-08-25
 ### Added
 - Added `State::spinodal` that calculates both spinodal points for a given temperature and composition using the same objective function that is also used in the critical point calculation. [#23](https://github.com/feos-org/feos/pull/23)

diff --git a/feos-core/src/parameter/mod.rs b/feos-core/src/parameter/mod.rs
@@ -241,7 +241,7 @@ where
         let n = pure_records.len();
         let mut binary_records = Array2::default([n, n]);
         for i in 0..n {
-            for j in 0..n {
+            for j in i + 1..n {
                 let mut vec = Vec::new();
                 for (id1, &n1) in segment_counts[i].iter() {
                     for (id2, &n2) in segment_counts[j].iter() {
@@ -253,7 +253,9 @@ where
                         vec.push((binary, n1, n2));
                     }
                 }
-                binary_records[(i, j)] = Self::Binary::from_segments_binary(&vec)?
+                let kij = Self::Binary::from_segments_binary(&vec)?;
+                binary_records[(i, j)] = kij.clone();
+                binary_records[(j, i)] = kij;
             }
         }
 

diff --git a/feos-core/src/parameter/model_record.rs b/feos-core/src/parameter/model_record.rs
@@ -2,7 +2,8 @@ use super::identifier::Identifier;
 use super::segment::SegmentRecord;
 use super::ParameterError;
 use conv::ValueInto;
-use serde::{Deserialize, Serialize};
+use serde::{de::DeserializeOwned, Deserialize, Serialize};
+use std::{fs::File, io::BufReader, path::Path};
 
 /// A collection of parameters of a pure substance.
 #[derive(Serialize, Deserialize, Debug, Clone)]
@@ -121,6 +122,15 @@ impl<I, B> BinaryRecord<I, B> {
             model_record,
         }
     }
+
+    /// Read a list of `BinaryRecord`s from a JSON file.
+    pub fn from_json<P: AsRef<Path>>(file: P) -> Result<Vec<Self>, ParameterError>
+    where
+        I: DeserializeOwned,
+        B: DeserializeOwned,
+    {
+        Ok(serde_json::from_reader(BufReader::new(File::open(file)?))?)
+    }
 }
 
 impl<I, B> std::fmt::Display for BinaryRecord<I, B>

diff --git a/feos-core/src/python/cubic.rs b/feos-core/src/python/cubic.rs
@@ -7,7 +7,7 @@ use crate::python::joback::PyJobackRecord;
 use crate::python::parameter::PyIdentifier;
 use crate::*;
 use ndarray::Array2;
-use numpy::PyReadonlyArray2;
+use numpy::{PyArray2, PyReadonlyArray2, ToPyArray};
 use pyo3::exceptions::PyTypeError;
 use pyo3::prelude::*;
 use std::convert::{TryFrom, TryInto};

diff --git a/feos-core/src/python/parameter.rs b/feos-core/src/python/parameter.rs
@@ -216,6 +216,25 @@ macro_rules! impl_binary_record {
                 }
             }
 
+            /// Read a list of `BinaryRecord`s from a JSON file.
+            ///
+            /// Parameters
+            /// ----------
+            /// path : str
+            ///     Path to file containing the binary records.
+            ///
+            /// Returns
+            /// -------
+            /// [BinaryRecord]
+            #[staticmethod]
+            #[pyo3(text_signature = "(path)")]
+            fn from_json(path: &str) -> Result<Vec<Self>, ParameterError> {
+                Ok(BinaryRecord::from_json(path)?
+                    .into_iter()
+                    .map(Self)
+                    .collect())
+            }
+
             #[getter]
             fn get_id1(&self) -> PyIdentifier {
                 PyIdentifier(self.0.id1.clone())
@@ -293,6 +312,25 @@ impl PyBinarySegmentRecord {
         Ok(Self(BinaryRecord::new(id1, id2, model_record)))
     }
 
+    /// Read a list of `BinarySegmentRecord`s from a JSON file.
+    ///
+    /// Parameters
+    /// ----------
+    /// path : str
+    ///     Path to file containing the binary records.
+    ///
+    /// Returns
+    /// -------
+    /// [BinarySegmentRecord]
+    #[staticmethod]
+    #[pyo3(text_signature = "(path)")]
+    fn from_json(path: &str) -> Result<Vec<Self>, ParameterError> {
+        Ok(BinaryRecord::from_json(path)?
+            .into_iter()
+            .map(Self)
+            .collect())
+    }
+
     #[getter]
     fn get_id1(&self) -> String {
         self.0.id1.clone()
@@ -470,12 +508,12 @@ macro_rules! impl_segment_record {
             ///
             /// Returns
             /// -------
-            /// SegmentRecord
+            /// [SegmentRecord]
             #[staticmethod]
             fn from_json(path: &str) -> Result<Vec<Self>, ParameterError> {
                 Ok(SegmentRecord::from_json(path)?
                     .into_iter()
-                    .map(|s| Self(s))
+                    .map(Self)
                     .collect())
             }
 
@@ -694,6 +732,11 @@ macro_rules! impl_parameter {
                     .map(|r| PyPureRecord(r.clone()))
                     .collect()
             }
+
+            #[getter]
+            fn get_binary_records<'py>(&self, py: Python<'py>) -> &'py PyArray2<f64> {
+                self.0.records().1.mapv(f64::from).view().to_pyarray(py)
+            }
         }
     };
 }

diff --git a/src/gc_pcsaft/eos/polar.rs b/src/gc_pcsaft/eos/polar.rs
@@ -68,19 +68,12 @@ impl Dipole {
         let ndipole = parameters.dipole_comp.len();
 
         let f2_term = Array2::from_shape_fn([ndipole; 2], |(i, j)| {
-            let di = parameters.dipole_comp[i];
-            let dj = parameters.dipole_comp[j];
-            parameters.mu2[i] * parameters.mu2[j] / parameters.s_ij[[di, dj]].powi(3)
+            parameters.mu2[i] * parameters.mu2[j] / parameters.s_ij[[i, j]].powi(3)
         });
 
         let f3_term = Array3::from_shape_fn([ndipole; 3], |(i, j, k)| {
-            let di = parameters.dipole_comp[i];
-            let dj = parameters.dipole_comp[j];
-            let dk = parameters.dipole_comp[k];
             parameters.mu2[i] * parameters.mu2[j] * parameters.mu2[k]
-                / (parameters.s_ij[[di, dj]]
-                    * parameters.s_ij[[di, dk]]
-                    * parameters.s_ij[[dj, dk]])
+                / (parameters.s_ij[[i, j]] * parameters.s_ij[[i, k]] * parameters.s_ij[[j, k]])
         });
 
         let mut mij1 = Array2::zeros((ndipole, ndipole));
@@ -131,38 +124,23 @@ impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for Dipole {
         let ndipole = p.dipole_comp.len();
 
         let t_inv = state.temperature.inv();
-        let eps_ij_t = Array2::from_shape_fn([ndipole; 2], |(i, j)| {
-            let di = p.dipole_comp[i];
-            let dj = p.dipole_comp[j];
-            t_inv * p.e_k_ij[[di, dj]]
-        });
+        let eps_ij_t = Array2::from_shape_fn([ndipole; 2], |(i, j)| t_inv * p.e_k_ij[[i, j]]);
 
         let rho = &state.partial_density;
-        let r = p.hs_diameter(state.temperature) * 0.5;
-        let eta = (rho * &p.m * &r * &r * &r).sum() * 4.0 * FRAC_PI_3;
+        let eta = p.zeta(state.temperature, &state.partial_density, [3])[0];
 
         let mut phi2 = D::zero();
         let mut phi3 = D::zero();
         for i in 0..ndipole {
             let di = p.dipole_comp[i];
             phi2 -= (rho[di] * rho[di] * self.f2_term[[i, i]])
-                * pair_integral_ij(
-                    self.mij1[[i, i]],
-                    self.mij2[[i, i]],
-                    eta,
-                    eps_ij_t[[di, di]],
-                );
+                * pair_integral_ij(self.mij1[[i, i]], self.mij2[[i, i]], eta, eps_ij_t[[i, i]]);
             phi3 -= (rho[di] * rho[di] * rho[di] * self.f3_term[[i, i, i]])
                 * triplet_integral_ijk(self.mijk1[[i, i, i]], self.mijk2[[i, i, i]], eta);
             for j in (i + 1)..ndipole {
                 let dj = p.dipole_comp[j];
                 phi2 -= (rho[di] * rho[dj] * self.f2_term[[i, j]])
-                    * pair_integral_ij(
-                        self.mij1[[i, j]],
-                        self.mij2[[i, j]],
-                        eta,
-                        eps_ij_t[[di, dj]],
-                    )
+                    * pair_integral_ij(self.mij1[[i, j]], self.mij2[[i, j]], eta, eps_ij_t[[i, j]])
                     * 2.0;
                 phi3 -= (rho[di] * rho[di] * rho[dj] * self.f3_term[[i, i, j]])
                     * triplet_integral_ijk(self.mijk1[[i, i, j]], self.mijk2[[i, i, j]], eta)

diff --git a/tests/gc_pcsaft/binary.rs b/tests/gc_pcsaft/binary.rs
@@ -3,9 +3,9 @@ use feos::gc_pcsaft::{GcPcSaft, GcPcSaftEosParameters};
 #[cfg(feature = "dft")]
 use feos::gc_pcsaft::{GcPcSaftFunctional, GcPcSaftFunctionalParameters};
 use feos_core::parameter::{IdentifierOption, ParameterHetero};
-use feos_core::{EosResult, State};
+use feos_core::{Contributions, EosResult, State};
 use ndarray::arr1;
-use quantity::si::{KELVIN, MOL};
+use quantity::si::{KELVIN, METER, MOL};
 use std::sync::Arc;
 
 #[test]
@@ -50,3 +50,31 @@ fn test_binary() -> EosResult<()> {
     );
     Ok(())
 }
+
+#[test]
+fn test_polar_term() -> EosResult<()> {
+    let parameters1 = GcPcSaftEosParameters::from_json_segments(
+        &["CCCOC(C)=O", "CCCO"],
+        "parameters/pcsaft/gc_substances.json",
+        "parameters/pcsaft/sauer2014_hetero.json",
+        None,
+        IdentifierOption::Smiles,
+    )?;
+    let parameters2 = GcPcSaftEosParameters::from_json_segments(
+        &["CCCO", "CCCOC(C)=O"],
+        "parameters/pcsaft/gc_substances.json",
+        "parameters/pcsaft/sauer2014_hetero.json",
+        None,
+        IdentifierOption::Smiles,
+    )?;
+    let eos1 = Arc::new(GcPcSaft::new(Arc::new(parameters1)));
+    let eos2 = Arc::new(GcPcSaft::new(Arc::new(parameters2)));
+    let moles = arr1(&[0.5, 0.5]) * MOL;
+    let p1 = State::new_nvt(&eos1, 300.0 * KELVIN, METER.powi(3), &moles)?
+        .pressure(Contributions::Total);
+    let p2 = State::new_nvt(&eos2, 300.0 * KELVIN, METER.powi(3), &moles)?
+        .pressure(Contributions::Total);
+    println!("{p1} {p2}");
+    assert_eq!(p1, p2);
+    Ok(())
+}