use crate::saftvrqmie::eos::FeynmanHibbsOrder;
use core::cmp::max;
use feos_core::parameter::Parameters;
use feos_core::{FeosError, FeosResult};
use nalgebra::{DMatrix, DVector};
use quantity::{KILOGRAM, NAV};
use serde::{Deserialize, Serialize};
use std::convert::TryFrom;
#[derive(Serialize, Deserialize, Clone, Default)]
pub struct SaftVRQMieRecord {
pub m: f64,
pub sigma: f64,
pub epsilon_k: f64,
pub lr: f64,
pub la: f64,
pub fh: usize,
}
impl SaftVRQMieRecord {
pub fn new(
m: f64,
sigma: f64,
epsilon_k: f64,
lr: f64,
la: f64,
fh: usize,
) -> FeosResult<SaftVRQMieRecord> {
if m != 1.0 {
return Err(FeosError::IncompatibleParameters(
"Segment number `m` is not one. Chain-contributions are currently not supported."
.to_string(),
));
}
Ok(SaftVRQMieRecord {
m,
sigma,
epsilon_k,
lr,
la,
fh,
})
}
}
#[derive(Serialize, Deserialize, Clone, Copy, Default)]
pub struct SaftVRQMieBinaryRecord {
pub k_ij: f64,
pub l_ij: f64,
}
pub type SaftVRQMieParameters = Parameters<SaftVRQMieRecord, SaftVRQMieBinaryRecord, ()>;
pub struct SaftVRQMiePars {
pub m: DVector<f64>,
pub sigma: DVector<f64>,
pub epsilon_k: DVector<f64>,
pub sigma_ij: DMatrix<f64>,
pub epsilon_k_ij: DMatrix<f64>,
pub c_ij: DMatrix<f64>,
pub lambda_r_ij: DMatrix<f64>,
pub lambda_a_ij: DMatrix<f64>,
pub mass_ij: DMatrix<f64>,
pub fh_ij: DMatrix<FeynmanHibbsOrder>,
}
impl SaftVRQMiePars {
pub fn new(parameters: &SaftVRQMieParameters) -> FeosResult<Self> {
let n = parameters.pure.len();
let [fh] = parameters.collate(|pr| [pr.fh]);
let [m, sigma, epsilon_k] = parameters.collate(|pr| [pr.m, pr.sigma, pr.epsilon_k]);
let [lr, la] = parameters.collate(|pr| [pr.lr, pr.la]);
let molarweight = ¶meters.molar_weight;
for (i, m) in m.iter().enumerate() {
if *m != 1.0 {
return Err(FeosError::IncompatibleParameters(format!(
"Segment number `m` for component {i} is not one. Chain-contributions are currently not supported."
)));
}
}
let mut fh_ij: DMatrix<FeynmanHibbsOrder> =
DMatrix::from_element(n, n, FeynmanHibbsOrder::FH0);
let [k_ij, l_ij] = parameters.collate_binary(|b| [b.k_ij, b.l_ij]);
let mut epsilon_k_ij = DMatrix::zeros(n, n);
let mut sigma_ij = DMatrix::zeros(n, n);
let mut e_k_ij = DMatrix::zeros(n, n);
let mut lambda_r_ij = DMatrix::zeros(n, n);
let mut lambda_a_ij = DMatrix::zeros(n, n);
let mut c_ij = DMatrix::zeros(n, n);
let mut mass_ij = DMatrix::zeros(n, n);
for i in 0..n {
for j in 0..n {
sigma_ij[(i, j)] = (1.0 - l_ij[(i, j)]) * 0.5 * (sigma[i] + sigma[j]);
e_k_ij[(i, j)] = (sigma[i].powi(3) * sigma[j].powi(3)).sqrt()
/ sigma_ij[(i, j)].powi(3)
* (epsilon_k[i] * epsilon_k[j]).sqrt();
epsilon_k_ij[(i, j)] = (1.0 - k_ij[(i, j)]) * e_k_ij[(i, j)];
lambda_r_ij[(i, j)] = ((lr[i] - 3.0) * (lr[j] - 3.0)).sqrt() + 3.0;
lambda_a_ij[(i, j)] = ((la[i] - 3.0) * (la[j] - 3.0)).sqrt() + 3.0;
c_ij[(i, j)] = lambda_r_ij[(i, j)] / (lambda_r_ij[(i, j)] - lambda_a_ij[(i, j)])
* (lambda_r_ij[(i, j)] / lambda_a_ij[(i, j)])
.powf(lambda_a_ij[(i, j)] / (lambda_r_ij[(i, j)] - lambda_a_ij[(i, j)]));
mass_ij[(i, j)] = (2.0 * molarweight.get(i) * molarweight.get(j)
/ (molarweight.get(i) + molarweight.get(j)))
.convert_into(KILOGRAM * NAV);
fh_ij[(i, j)] = FeynmanHibbsOrder::try_from(max(fh[i], fh[j]))?;
if fh[i] * fh[j] == 2 {
return Err(FeosError::IncompatibleParameters(format!(
"cannot combine Feynman-Hibbs orders 1 and 2. Component {} has order {} and component {} has order {}.",
i, fh[i], j, fh[j]
)));
}
}
}
Ok(Self {
m,
sigma,
epsilon_k,
sigma_ij,
epsilon_k_ij,
c_ij,
lambda_r_ij,
lambda_a_ij,
mass_ij,
fh_ij,
})
}
}
#[cfg(test)]
pub mod utils {
use feos_core::parameter::PureRecord;
use super::*;
pub fn hydrogen_fh(fh: &str) -> SaftVRQMiePars {
let hydrogen_json = &(r#"
{
"identifier": {
"cas": "1333-74-0",
"name": "hydrogen",
"iupac_name": "hydrogen",
"smiles": "[HH]",
"inchi": "InChI=1S/H2/h1H",
"formula": "H2"
},
"m": 1.0,
"sigma": 3.0243,
"epsilon_k": 26.706,
"lr": 9.0,
"la": 6.0,
"fh": "#
.to_owned()
+ fh
+ r#",
"molarweight": 2.0157309551872
}"#);
let hydrogen_record: PureRecord<SaftVRQMieRecord, ()> =
serde_json::from_str(hydrogen_json).expect("Unable to parse json.");
SaftVRQMiePars::new(&SaftVRQMieParameters::new_pure(hydrogen_record).unwrap()).unwrap()
}
pub fn helium_fh1() -> PureRecord<SaftVRQMieRecord, ()> {
let helium_json = r#"
{
"identifier": {
"cas": "1333-74-0",
"name": "helium",
"iupac_name": "helium",
"smiles": "[HH]",
"inchi": "InChI=1S/H2/h1H",
"formula": "He"
},
"m": 1.0,
"sigma": 2.7443,
"epsilon_k": 5.4195,
"lr": 9.0,
"la": 6.0,
"fh": 1,
"molarweight": 4.002601643881807
}"#;
serde_json::from_str(helium_json).expect("Unable to parse json.")
}
pub fn hydrogen_fh2() -> PureRecord<SaftVRQMieRecord, ()> {
let helium_json = r#"
{
"identifier": {
"cas": "1333-74-0",
"name": "hydrogen",
"iupac_name": "hydrogen",
"smiles": "[HH]",
"inchi": "InChI=1S/H2/h1H",
"formula": "H2"
},
"m": 1.0,
"sigma": 3.0243,
"epsilon_k": 26.706,
"lr": 9.0,
"la": 6.0,
"fh": 2,
"molarweight": 2.0157309551872
}"#;
serde_json::from_str(helium_json).expect("Unable to parse json.")
}
#[expect(dead_code)]
pub fn neon_fh1() -> SaftVRQMieParameters {
let neon_json = r#"
{
"identifier": {
"cas": "1333-74-0",
"name": "neon",
"iupac_name": "neon",
"smiles": "[HH]",
"inchi": "InChI=1S/H2/h1H",
"formula": "Ne"
},
"m": 1.0,
"sigma": 2.7778,
"epsilon_k": 37.501,
"lr": 13.0,
"la": 6.0,
"fh": 1,
"molarweight": 20.17969806457545
}"#;
let neon_record: PureRecord<SaftVRQMieRecord, ()> =
serde_json::from_str(neon_json).expect("Unable to parse json.");
SaftVRQMieParameters::new_pure(neon_record).unwrap()
}
pub fn h2_ne_fh(fh: &str) -> SaftVRQMiePars {
let binary_json = &(r#"[
{
"identifier": {
"cas": "1333-74-0",
"name": "hydrogen",
"iupac_name": "hydrogen",
"smiles": "[HH]",
"inchi": "InChI=1S/H2/h1H",
"formula": "H2"
},
"m": 1.0,
"sigma": 3.0243,
"epsilon_k": 26.706,
"lr": 9.0,
"la": 6.0,
"fh": "#
.to_owned()
+ fh
+ r#",
"molarweight": 2.0157309551872
},
{
"identifier": {
"cas": "1333-74-0",
"name": "neon",
"iupac_name": "neon",
"smiles": "[HH]",
"inchi": "InChI=1S/H2/h1H",
"formula": "H2"
},
"m": 1.0,
"sigma": 2.7778,
"epsilon_k": 37.501,
"lr": 13.0,
"la": 6.0,
"fh": "#
+ fh
+ r#",
"molarweight": 20.17969806457545
}
]"#);
let binary_record: [PureRecord<SaftVRQMieRecord, ()>; 2] =
serde_json::from_str(binary_json).expect("Unable to parse json.");
SaftVRQMiePars::new(
&SaftVRQMieParameters::new_binary(
binary_record,
Some(SaftVRQMieBinaryRecord {
k_ij: 0.105,
l_ij: 0.0,
}),
vec![],
)
.unwrap(),
)
.unwrap()
}
}
#[cfg(test)]
mod test {
use super::SaftVRQMieParameters;
use super::utils::{helium_fh1, hydrogen_fh2};
use crate::saftvrqmie::SaftVRQMie;
#[test]
#[should_panic(
expected = "cannot combine Feynman-Hibbs orders 1 and 2. Component 0 has order 1 and component 1 has order 2."
)]
fn incompatible_order() {
let order1 = helium_fh1();
let order2 = hydrogen_fh2();
SaftVRQMie::new(SaftVRQMieParameters::new_binary([order1, order2], None, vec![]).unwrap())
.unwrap();
}
}