From bf56d5c328ad686f9581396ead0fd2b7c1ab070e Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 7 Jan 2023 13:41:02 +0100 Subject: [PATCH 1/5] Group contribution method for binary interaction parameters --- feos-core/CHANGELOG.md | 7 +++- feos-core/src/parameter/mod.rs | 6 ++-- feos-core/src/parameter/model_record.rs | 12 ++++++- feos-core/src/python/cubic.rs | 2 +- feos-core/src/python/parameter.rs | 47 +++++++++++++++++++++++-- src/gc_pcsaft/eos/polar.rs | 34 ++++-------------- tests/gc_pcsaft/binary.rs | 32 +++++++++++++++-- 7 files changed, 103 insertions(+), 37 deletions(-) diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md index 55d11ddd1..3137cd4c1 100644 --- 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 index f592b214f..48c62c303 100644 --- 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 index baa34a8a8..876ac0589 100644 --- 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 BinaryRecord { model_record, } } + + /// Read a list of `BinaryRecord`s from a JSON file. + pub fn from_json>(file: P) -> Result, ParameterError> + where + I: DeserializeOwned, + B: DeserializeOwned, + { + Ok(serde_json::from_reader(BufReader::new(File::open(file)?))?) + } } impl std::fmt::Display for BinaryRecord diff --git a/feos-core/src/python/cubic.rs b/feos-core/src/python/cubic.rs index 8096842d3..80a82b3f7 100644 --- 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 index adb90ce69..a6dda5884 100644 --- 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, 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, 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, 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 { + 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 index 9256524d1..e72dc2f37 100644 --- 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> HelmholtzEnergyDual 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 index a9430ca7c..d3a1dab45 100644 --- 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(()) +} From 89063d76da07be8de8da5d3452563fbb1d4c9dfa Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 7 Jan 2023 13:42:46 +0100 Subject: [PATCH 2/5] add changelog entry --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 22cfe19f2..c81da85d5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated `num-dual` dependency to 0.6. [#119](https://github.com/feos-org/feos/pull/119) +### Fixed +- Fixed incorrect indexing that lead to panics in the polar contribution of gc-PC-SAFT. [#104](https://github.com/feos-org/feos/pull/104) + ## [0.3.0] - 2022-09-14 - Major restructuring of the entire `feos` project. All individual models are reunited in the `feos` crate. `feos-core` and `feos-dft` still live as individual crates within the `feos` workspace. From f20a1ec9b19c02c3ebe774a6b8322239e5a6d14b Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 7 Jan 2023 13:44:59 +0100 Subject: [PATCH 3/5] not a big fan of how rust-analyzer does imports --- feos-core/src/parameter/model_record.rs | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/feos-core/src/parameter/model_record.rs b/feos-core/src/parameter/model_record.rs index 876ac0589..382b68936 100644 --- a/feos-core/src/parameter/model_record.rs +++ b/feos-core/src/parameter/model_record.rs @@ -2,8 +2,11 @@ use super::identifier::Identifier; use super::segment::SegmentRecord; use super::ParameterError; use conv::ValueInto; -use serde::{de::DeserializeOwned, Deserialize, Serialize}; -use std::{fs::File, io::BufReader, path::Path}; +use serde::de::DeserializeOwned; +use serde::{Deserialize, Serialize}; +use std::fs::File; +use std::io::BufReader; +use std::path::Path; /// A collection of parameters of a pure substance. #[derive(Serialize, Deserialize, Debug, Clone)] From c056a0123caa2a16304d5a3958c191c9cc7d61d0 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 7 Jan 2023 13:57:28 +0100 Subject: [PATCH 4/5] Fixes for UV and SAFTVRQMie --- feos-core/src/python/parameter.rs | 7 ++++++- src/uvtheory/python.rs | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/feos-core/src/python/parameter.rs b/feos-core/src/python/parameter.rs index a6dda5884..50fb8592a 100644 --- a/feos-core/src/python/parameter.rs +++ b/feos-core/src/python/parameter.rs @@ -735,7 +735,12 @@ macro_rules! impl_parameter { #[getter] fn get_binary_records<'py>(&self, py: Python<'py>) -> &'py PyArray2 { - self.0.records().1.mapv(f64::from).view().to_pyarray(py) + self.0 + .records() + .1 + .mapv(|r| f64::try_from(r).unwrap()) + .view() + .to_pyarray(py) } } }; diff --git a/src/uvtheory/python.rs b/src/uvtheory/python.rs index 08c135a18..fa235b6a5 100644 --- a/src/uvtheory/python.rs +++ b/src/uvtheory/python.rs @@ -6,7 +6,7 @@ use feos_core::parameter::{ use feos_core::python::parameter::*; use feos_core::*; use ndarray::Array2; -use numpy::PyReadonlyArray2; +use numpy::{PyArray2, PyReadonlyArray2, ToPyArray}; use pyo3::exceptions::PyTypeError; use pyo3::prelude::*; use std::convert::{TryFrom, TryInto}; From 01348c7a771ce09e730081691c58f250a95adf51 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 27 Jan 2023 09:44:44 +0100 Subject: [PATCH 5/5] other combination rule and test for kij --- src/pcsaft/parameters.rs | 45 +++++++++++++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 9517b9317..7053aa6a6 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -255,11 +255,11 @@ impl From for f64 { impl> FromSegmentsBinary for PcSaftBinaryRecord { fn from_segments_binary(segments: &[(Self, T, T)]) -> Result { - let k_ij = segments - .iter() - .map(|(br, n1, n2)| br.k_ij * (*n1).value_into().unwrap() * (*n2).value_into().unwrap()) - .sum(); - Ok(Self { k_ij }) + let (k_ij, n) = segments.iter().fold((0.0, 0.0), |(k_ij, n), (br, n1, n2)| { + let nab = (*n1).value_into().unwrap() * (*n2).value_into().unwrap(); + (k_ij + br.k_ij * nab, n + nab) + }); + Ok(Self { k_ij: k_ij / n }) } } @@ -523,6 +523,7 @@ impl std::fmt::Display for PcSaftParameters { pub mod utils { use super::*; use feos_core::joback::JobackRecord; + use feos_core::parameter::{BinaryRecord, ChemicalRecord, SegmentRecord}; use std::sync::Arc; pub fn propane_parameters() -> Arc { @@ -733,4 +734,38 @@ pub mod utils { serde_json::from_str(binary_json).expect("Unable to parse json."); Arc::new(PcSaftParameters::new_binary(binary_record, None)) } + + #[test] + pub fn test_kij() -> Result<(), ParameterError> { + let ch3: String = "CH3".into(); + let ch2: String = "CH2".into(); + let oh = "OH".into(); + let propane = ChemicalRecord::new( + Default::default(), + vec![ch3.clone(), ch2.clone(), ch3.clone()], + None, + ); + let ethanol = ChemicalRecord::new(Default::default(), vec![ch3, ch2, oh], None); + let segment_records = SegmentRecord::from_json("parameters/pcsaft/sauer2014_homo.json")?; + let kij = [("CH3", "OH", -0.2), ("CH2", "OH", -0.1)]; + let binary_segment_records = kij + .iter() + .map(|&(id1, id2, k_ij)| { + BinaryRecord::new(id1.into(), id2.into(), PcSaftBinaryRecord { k_ij }) + }) + .collect(); + let params = PcSaftParameters::from_segments( + vec![propane, ethanol], + segment_records, + Some(binary_segment_records), + )?; + let k_ij = ¶ms.binary_records; + assert_eq!(k_ij[[0, 0]].k_ij, 0.0); + assert_eq!(k_ij[[0, 1]].k_ij, -0.5 / 9.); + assert_eq!(k_ij[[1, 0]].k_ij, -0.5 / 9.); + assert_eq!(k_ij[[1, 1]].k_ij, 0.0); + println!("{k_ij}"); + + Ok(()) + } }