feos_ad/eos/
gc_pcsaft.rs

1use super::pcsaft::{A0, A1, A2, AD, B0, B1, B2, BD, CD};
2use crate::{ParametersAD, ResidualHelmholtzEnergy};
3use nalgebra::{SMatrix, SVector};
4use num_dual::DualNum;
5use quantity::{JOULE, KB, KELVIN};
6use std::collections::HashMap;
7use std::f64::consts::{FRAC_PI_6, PI};
8use std::sync::LazyLock;
9
10const PI_SQ_43: f64 = 4.0 / 3.0 * PI * PI;
11
12const MAX_ETA: f64 = 0.5;
13
14const N_GROUPS: usize = 20;
15const GROUPS: [&str; N_GROUPS] = [
16    "CH3", "CH2", ">CH", ">C<", "=CH2", "=CH", "=C<", "C≡CH", "CH2_hex", "CH_hex", "CH2_pent",
17    "CH_pent", "CH_arom", "C_arom", "CH=O", ">C=O", "OCH3", "OCH2", "HCOO", "COO",
18];
19const M: [f64; N_GROUPS] = [
20    0.77247, 0.7912, 0.52235, -0.70131, 0.70581, 0.90182, 0.98505, 1.1615, 0.8793, 0.42115,
21    0.90057, 0.69343, 0.88259, 0.77531, 1.1889, 1.1889, 1.1907, 1.1817, 1.2789, 1.2869,
22];
23const SIGMA: [f64; N_GROUPS] = [
24    3.6937, 3.0207, 0.99912, 0.5435, 3.163, 2.8864, 2.245, 3.3187, 2.9995, 1.3078, 3.0437, 1.2894,
25    2.9475, 1.6719, 3.2948, 3.1026, 2.7795, 3.009, 3.373, 3.0643,
26];
27const EPSILON_K: [f64; N_GROUPS] = [
28    181.49, 157.23, 269.84, 0.0, 171.34, 158.9, 146.86, 255.13, 157.93, 131.79, 158.34, 140.69,
29    156.51, 178.81, 316.91, 280.43, 284.91, 203.11, 307.44, 273.9,
30];
31const MU: [f64; N_GROUPS] = [
32    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.4126, 3.4167, 0.0,
33    2.6945, 2.6808, 3.3428,
34];
35
36const K_AB_LIST: [((&str, &str), f64); 128] = [
37    (("CH3", "C≡CH"), 0.1652754031504817),
38    (("CH2", "C≡CH"), -0.006553805247652898),
39    ((">CH", "C≡CH"), -0.45179562711827204),
40    ((">CH", "CH_arom"), -0.2623133986577006),
41    (("CH3", "C_arom"), 0.1029438146806397),
42    (("CH2", "C_arom"), -0.019417887747868838),
43    (("CH3", "CH_arom"), 0.03524440115729556),
44    (("CH2", "CH_arom"), -0.005376660474268035),
45    ((">CH", "C_arom"), -0.4536159179058391),
46    (("CH3", "CH=O"), 0.07114002498453002),
47    (("CH2", "CH=O"), 0.05810424929073859),
48    ((">C=O", "CH3"), 0.06638452980811302),
49    ((">C=O", "CH2"), 0.053964852486994674),
50    ((">C=O", ">CH"), -0.01886821792425231),
51    (("CH3", "OCH2"), 0.0040270442665658705),
52    ((">CH", "OCH3"), -0.217389127174595),
53    (("CH2", "OCH2"), -0.015001007000489858),
54    (("CH3", "OCH3"), 0.04370031373865),
55    ((">CH", "OCH2"), -0.038601031822155345),
56    (("CH2", "OCH3"), 0.021857618541656035),
57    (("CH3", "HCOO"), 0.046221094393257035),
58    (("CH2", "HCOO"), 0.07467469076976556),
59    ((">CH", "COO"), 0.17498779403987508),
60    (("CH3", "COO"), 0.05050525601666908),
61    (("CH2", "COO"), 0.04126150329459275),
62    (("=CH2", "C≡CH"), 0.28518205265981117),
63    (("=CH", "C≡CH"), -0.19471892756822953),
64    (("=C<", "C≡CH"), -0.37966037059629637),
65    (("=CH", "C_arom"), -0.04794957652572591),
66    (("=C<", "C_arom"), -0.19059760666493372),
67    (("=CH2", "C_arom"), 0.002829235390387696),
68    (("=CH", "CH_arom"), 0.008744506522926979),
69    (("=C<", "CH_arom"), -0.15988275905975652),
70    (("=CH2", "CH_arom"), 0.039714702972518216),
71    (("=C<", ">C=O"), -0.23303936628601363),
72    (("=CH2", ">C=O"), 0.0451103450863995),
73    (("=CH", ">C=O"), 0.0028282796817093118),
74    (("=CH2", "OCH3"), -0.017336555493171858),
75    (("=CH", "OCH2"), 0.028713611730255537),
76    (("=C<", "OCH2"), -0.040576835006969125),
77    (("=CH2", "OCH2"), 0.02792348761278379),
78    (("=CH", "OCH3"), 0.021854680107075346),
79    (("=C<", "OCH3"), -0.21464227012985213),
80    (("=CH", "HCOO"), -0.021573291908933357),
81    (("=C<", "HCOO"), -0.021791386613505864),
82    (("=CH2", "COO"), -0.08063693709595326),
83    (("=CH", "COO"), -0.07829355920744586),
84    (("=C<", "COO"), -0.19510136763283895),
85    (("CH_arom", "C≡CH"), -0.04955767628386867),
86    (("C_arom", "C≡CH"), -0.04953394596854589),
87    (("CH=O", "C≡CH"), -0.33948211818518437),
88    ((">C=O", "C≡CH"), -0.3657376137845608),
89    (("C≡CH", "OCH2"), -0.3344648388007797),
90    (("C≡CH", "OCH3"), -0.38290586519600983),
91    (("C≡CH", "HCOO"), -0.24272079727170506),
92    (("COO", "C≡CH"), -0.3654438081227738),
93    (("C≡CH", "OH"), -0.10409652963367089),
94    (("CH_arom", "C_arom"), -0.12728867698722554),
95    (("CH_arom", "CH_arom"), -0.023433119170883764),
96    (("C_arom", "C_arom"), -0.28421918877688607),
97    (("CH2_hex", "CH_arom"), 0.001447550584366187),
98    (("CH_hex", "C_arom"), -0.40316115168074723),
99    (("CH_arom", "CH_hex"), -0.3576321797883022),
100    (("CH2_hex", "C_arom"), 0.08521616156959391),
101    (("CH_arom", "CH_pent"), -0.18714862430161244),
102    (("CH_pent", "C_arom"), 0.03457714936255159),
103    (("CH2_pent", "C_arom"), 0.03812116290565423),
104    (("CH2_pent", "CH_arom"), 0.004966663517893626),
105    (("CH=O", "C_arom"), -0.025044677100339502),
106    (("CH=O", "CH_arom"), -0.06309730593619743),
107    ((">C=O", "CH_arom"), -0.1106708081496441),
108    ((">C=O", "C_arom"), -0.04717854573415193),
109    (("C_arom", "OCH3"), -0.14006221692396298),
110    (("CH_arom", "OCH2"), -0.11744345395130776),
111    (("CH_arom", "OCH3"), -0.06563917699710337),
112    (("C_arom", "OCH2"), -0.03747169315781876),
113    (("CH_arom", "HCOO"), -0.04404325532489947),
114    (("C_arom", "HCOO"), 0.2898776740748664),
115    (("CH_arom", "COO"), -0.10600926342624745),
116    (("COO", "C_arom"), -0.11054573554364296),
117    (("C_arom", "OH"), 0.2990171120344822),
118    (("CH_arom", "OH"), -0.039398695604311314),
119    (("C_arom", "NH2"), 0.4535791221221567),
120    (("CH_arom", "NH2"), -0.06290638692043257),
121    (("CH2_hex", "CH=O"), 0.09047030071006708),
122    (("CH=O", "CH_hex"), -0.14747598417210014),
123    ((">C=O", "CH_hex"), 0.0676825668907787),
124    ((">C=O", "CH2_hex"), 0.09082375748804353),
125    (("CH2_hex", "OCH3"), 0.042823701076275526),
126    (("CH_hex", "OCH2"), -0.0936451919984422),
127    (("CH_hex", "OCH3"), 0.12111386208387202),
128    (("CH2_hex", "OCH2"), 0.013698887705260425),
129    (("CH2_hex", "HCOO"), 0.08719198819954514),
130    (("CH2_hex", "COO"), 0.05937938878778157),
131    (("CH_hex", "COO"), -0.10319900739370075),
132    (("CH2_hex", "OH"), 0.06127398203560399),
133    (("CH_hex", "OH"), 0.35825180807831797),
134    (("CH2_pent", "COO"), 0.05124310486288829),
135    (("CH_pent", "OH"), 0.11518421254437769),
136    (("CH2_pent", "OH"), 0.08215868571093943),
137    (("CH=O", "CH=O"), -0.1570556622280003),
138    ((">C=O", "CH=O"), -0.16452206370456918),
139    (("CH=O", "OCH2"), 0.0027095251191336708),
140    (("CH=O", "HCOO"), -0.07673278642721204),
141    (("CH=O", "COO"), -0.16365969940991995),
142    (("CH=O", "OH"), -0.12245349842770242),
143    ((">C=O", ">C=O"), -0.17931771307891634),
144    ((">C=O", "OCH3"), -0.20038719736411056),
145    ((">C=O", "OCH2"), 0.04468736703539099),
146    ((">C=O", "HCOO"), -0.13541978245760022),
147    ((">C=O", "COO"), -0.14605212162381323),
148    ((">C=O", "OH"), -0.1392769563372809),
149    ((">C=O", "NH2"), -0.371931948310995),
150    (("OCH2", "OCH2"), 0.09662571077941844),
151    (("OCH2", "OCH3"), -0.2812620283200189),
152    (("OCH3", "OCH3"), -0.13909723652059505),
153    (("HCOO", "OCH3"), -0.0929570619422749),
154    (("COO", "OCH2"), -0.11408007406963222),
155    (("COO", "OCH3"), -0.21710938244245623),
156    (("OCH2", "OH"), -0.014272196525878467),
157    (("OCH3", "OH"), -0.039585706351111166),
158    (("HCOO", "HCOO"), -0.14303475853269773),
159    (("COO", "HCOO"), -0.14056680898820434),
160    (("HCOO", "OH"), -0.11204049889700908),
161    (("COO", "COO"), -0.1879219131496382),
162    (("COO", "OH"), -0.09507071103459414),
163    (("COO", "NH2"), -0.2799573216348791),
164    (("NH2", "OH"), -0.42107448986356144),
165];
166
167static K_AB_MAP: LazyLock<HashMap<(&str, &str), f64>> = LazyLock::new(|| {
168    K_AB_LIST
169        .into_iter()
170        .map(|((s1, s2), val)| ((s2, s1), val))
171        .chain(K_AB_LIST)
172        .collect()
173});
174
175static K_AB: LazyLock<SMatrix<f64, N_GROUPS, N_GROUPS>> = LazyLock::new(|| {
176    SMatrix::from_fn(|i, j| *K_AB_MAP.get(&(GROUPS[i], GROUPS[j])).unwrap_or(&0.0))
177});
178
179/// Parameters used to instantiate [GcPcSaft].
180#[derive(Clone)]
181pub struct GcPcSaftParameters<D, const N: usize> {
182    pub groups: SMatrix<D, N_GROUPS, N>,
183    pub bonds: [Vec<([usize; 2], D)>; N],
184}
185
186impl<D: DualNum<f64> + Copy, const N: usize> GcPcSaftParameters<D, N> {
187    pub fn re(&self) -> GcPcSaftParameters<f64, N> {
188        let Self { groups, bonds } = self;
189        let groups = groups.map(|g| g.re());
190        let bonds = bonds
191            .each_ref()
192            .map(|b| b.iter().map(|&(b, v)| (b, v.re())).collect());
193        GcPcSaftParameters { groups, bonds }
194    }
195}
196
197impl<D: DualNum<f64> + Copy, const N: usize> GcPcSaftParameters<D, N> {
198    pub fn from_groups(
199        group_map: [&HashMap<&'static str, D>; N],
200        bond_map: [&HashMap<[&'static str; 2], D>; N],
201    ) -> Self {
202        let groups =
203            SMatrix::from(group_map.map(|r| GROUPS.map(|g| *r.get(g).unwrap_or(&D::zero()))));
204        let group_indices: HashMap<_, _> = GROUPS
205            .into_iter()
206            .enumerate()
207            .map(|(g, i)| (i, g))
208            .collect();
209        let bonds = bond_map.map(|r| {
210            r.iter()
211                .map(|([g1, g2], &c)| ([group_indices[g1], group_indices[g2]], c))
212                .collect()
213        });
214        Self { groups, bonds }
215    }
216}
217
218/// The heterosegmented GC model for PC-SAFT by Sauer et al.
219pub struct GcPcSaft<const N: usize>(pub GcPcSaftParameters<f64, N>);
220
221impl<const N: usize> GcPcSaft<N> {
222    fn helmholtz_energy_density<D: DualNum<f64> + Copy>(
223        parameters: &GcPcSaftParameters<D, N>,
224        temperature: D,
225        density: &SVector<D, N>,
226    ) -> D {
227        let GcPcSaftParameters { groups, bonds } = parameters;
228
229        // convert parameters
230        let m = apply_group_count(groups, &SVector::from(M.map(D::from)));
231        let sigma = SVector::from(SIGMA.map(D::from));
232        let epsilon_k = SVector::from(EPSILON_K.map(D::from));
233        let mu = SVector::from(MU.map(D::from));
234
235        // temperature dependent segment diameter
236        let t_inv = temperature.recip();
237        let diameter = (epsilon_k * (t_inv * -3.0))
238            .map(|x| -x.exp() * 0.12 + 1.0)
239            .component_mul(&sigma);
240
241        // packing fractions
242        let mut zeta = [D::zero(); 4];
243        for c in 0..N {
244            for i in 0..diameter.len() {
245                for (z, &k) in zeta.iter_mut().zip([0, 1, 2, 3].iter()) {
246                    *z += density[c] * diameter[i].powi(k) * m[(i, c)] * FRAC_PI_6;
247                }
248            }
249        }
250        let zeta_23 = zeta[2] / zeta[3];
251        let frac_1mz3 = -(zeta[3] - 1.0).recip();
252
253        // hard sphere
254        let hs = (zeta[1] * zeta[2] * frac_1mz3 * 3.0
255            + zeta[2].powi(2) * frac_1mz3.powi(2) * zeta_23
256            + (zeta[2] * zeta_23.powi(2) - zeta[0]) * (zeta[3] * (-1.0)).ln_1p())
257            / std::f64::consts::FRAC_PI_6;
258
259        // hard chain
260        let c = zeta[2] * frac_1mz3 * frac_1mz3;
261        let hc: D = bonds
262            .iter()
263            .zip(density.iter())
264            .flat_map(|(bonds, &rho)| {
265                bonds.iter().map(move |([i, j], count)| {
266                    let (di, dj) = (diameter[*i], diameter[*j]);
267                    let cdij = c * di * dj / (di + dj);
268                    let g = frac_1mz3 + cdij * 3.0 - cdij * cdij * (zeta[3] - 1.0) * 2.0;
269                    -rho * count * g.ln()
270                })
271            })
272            .sum();
273
274        // packing fraction
275        let eta = zeta[3];
276
277        // mean segment number
278        let molefracs = density / density.sum();
279        let mbar = m.row_sum().tr_dot(&molefracs);
280
281        // crosswise interactions of all groups on all chains
282        let eps_ij = SVector::from(EPSILON_K).map(f64::sqrt);
283        let eps_ij = eps_ij * eps_ij.transpose();
284        let mut rho1mix = D::zero();
285        let mut rho2mix = D::zero();
286        for (m_i, &rho_i) in m.column_iter().zip(density.iter()) {
287            for (m_j, &rho_j) in m.column_iter().zip(density.iter()) {
288                for i in 0..N_GROUPS {
289                    for j in 0..N_GROUPS {
290                        let k_ab = if m_i != m_j { K_AB[(i, j)] } else { 0.0 };
291                        let eps_ij_t = temperature.recip() * eps_ij[(i, j)] * (1.0 - k_ab);
292                        let sigma_ij = ((SIGMA[i] + SIGMA[j]) * 0.5).powi(3);
293                        let rho1 = rho_i * rho_j * (eps_ij_t * m_i[i] * m_j[j] * sigma_ij);
294                        rho1mix += rho1;
295                        rho2mix += rho1 * eps_ij_t;
296                    }
297                }
298            }
299        }
300
301        // I1, I2 and C1
302        let mut i1 = D::zero();
303        let mut i2 = D::zero();
304        let mut eta_i = D::one();
305        let m1 = (mbar - 1.0) / mbar;
306        let m2 = (mbar - 2.0) / mbar * m1;
307        for i in 0..=6 {
308            i1 += (m2 * A2[i] + m1 * A1[i] + A0[i]) * eta_i;
309            i2 += (m2 * B2[i] + m1 * B1[i] + B0[i]) * eta_i;
310            eta_i *= eta;
311        }
312        let c1 = (mbar * (eta * 8.0 - eta.powi(2) * 2.0) / (eta - 1.0).powi(4)
313            + (D::one() - mbar)
314                * (eta * 20.0 - eta.powi(2) * 27.0 + eta.powi(3) * 12.0 - eta.powi(4) * 2.0)
315                / ((eta - 1.0) * (eta - 2.0)).powi(2)
316            + 1.0)
317            .recip();
318
319        // dispersion
320        let disp = (-rho1mix * i1 * 2.0 - rho2mix * mbar * c1 * i2) * PI;
321
322        // dipoles
323        let m_mix = m.row_sum();
324        let sigma_mix = apply_group_count(&m, &sigma.map(|s| s.powi(3)))
325            .row_sum()
326            .component_div(&m_mix)
327            .map(|s| s.cbrt());
328        let epsilon_k_mix = apply_group_count(&m, &epsilon_k)
329            .row_sum()
330            .component_div(&m_mix);
331        let mu2 = apply_group_count(groups, &mu.component_mul(&mu))
332            .row_sum()
333            .component_div(&m_mix)
334            * (t_inv * 1e-19 * JOULE.convert_into(KELVIN * KB));
335        let m_mix = m_mix.map(|m| if m.re() > 2.0 { D::from(2.0) } else { m });
336
337        let mut phi2 = D::zero();
338        let mut phi3 = D::zero();
339        for i in 0..N {
340            for j in i..N {
341                let sigma_ij_3 = ((sigma_mix[i] + sigma_mix[j]) * 0.5).powi(3);
342                let mij = (m_mix[i] * m_mix[j]).sqrt();
343                let mij1 = (mij - 1.0) / mij;
344                let mij2 = mij1 * (mij - 2.0) / mij;
345                let eps_ij_t = t_inv * (epsilon_k_mix[i] * epsilon_k_mix[j]).sqrt();
346                let c = if i == j { 1.0 } else { 2.0 };
347                phi2 -= (density[i] * density[j] * mu2[i] * mu2[j])
348                    * pair_integral(mij1, mij2, eta, eps_ij_t)
349                    / sigma_ij_3
350                    * c;
351                for k in j..N {
352                    let sigma_ij = (sigma_mix[i] + sigma_mix[j]) * 0.5;
353                    let sigma_ik = (sigma_mix[i] + sigma_mix[k]) * 0.5;
354                    let sigma_jk = (sigma_mix[j] + sigma_mix[k]) * 0.5;
355                    let mijk = (m_mix[i] * m_mix[j] * m_mix[k]).cbrt();
356                    let mijk1 = (mijk - 1.0) / mijk;
357                    let mijk2 = mijk1 * (mijk - 2.0) / mijk;
358                    let c = if i == j && i == k {
359                        1.0
360                    } else if i == j || i == k || j == k {
361                        3.0
362                    } else {
363                        6.0
364                    };
365                    phi3 -= (density[i] * density[j] * density[k] * mu2[i] * mu2[j] * mu2[k])
366                        * triplet_integral(mijk1, mijk2, eta)
367                        / (sigma_ij * sigma_ik * sigma_jk)
368                        * c;
369                }
370            }
371        }
372        phi2 *= PI;
373        phi3 *= PI_SQ_43;
374        let mut dipole = phi2 * phi2 / (phi2 - phi3);
375        if dipole.re().is_nan() {
376            dipole = phi2
377        }
378
379        (hs + hc + disp + dipole) * temperature
380    }
381}
382
383impl<const N: usize> ParametersAD for GcPcSaft<N> {
384    type Parameters<D: DualNum<f64> + Copy> = GcPcSaftParameters<D, N>;
385
386    fn params<D: DualNum<f64> + Copy>(&self) -> GcPcSaftParameters<D, N> {
387        let GcPcSaftParameters { groups, bonds } = &self.0;
388        let groups = groups.map(D::from);
389        let bonds = bonds
390            .each_ref()
391            .map(|b| b.iter().map(|&(b, v)| (b, D::from(v))).collect());
392        GcPcSaftParameters { groups, bonds }
393    }
394
395    fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy + Copy>(
396        parameters: &Self::Parameters<D>,
397    ) -> Self::Parameters<D2> {
398        let GcPcSaftParameters { groups, bonds } = parameters;
399        let groups = groups.map(D2::from_inner);
400        let bonds = bonds
401            .each_ref()
402            .map(|b| b.iter().map(|&(b, v)| (b, D2::from_inner(v))).collect());
403        GcPcSaftParameters { groups, bonds }
404    }
405}
406
407impl<const N: usize> ResidualHelmholtzEnergy<N> for GcPcSaft<N> {
408    const RESIDUAL: &str = "gc-PC-SAFT";
409
410    fn compute_max_density(&self, molefracs: &SVector<f64, N>) -> f64 {
411        let msigma3 = SVector::from_fn(|i, _| M[i] * SIGMA[i].powi(3));
412        let GcPcSaftParameters { groups, bonds: _ } = &self.0;
413        let msigma3 = apply_group_count(groups, &msigma3).row_sum();
414        let x: f64 = msigma3
415            .into_iter()
416            .zip(molefracs)
417            .map(|(ms3, x)| x * ms3)
418            .sum();
419        MAX_ETA / (x * FRAC_PI_6)
420    }
421
422    fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
423        parameters: &GcPcSaftParameters<D, N>,
424        temperature: D,
425        density: &SVector<D, N>,
426    ) -> D {
427        Self::helmholtz_energy_density(parameters, temperature, density)
428    }
429}
430
431fn apply_group_count<D: DualNum<f64> + Copy, const N: usize>(
432    groups: &SMatrix<D, N_GROUPS, N>,
433    x: &SVector<D, N_GROUPS>,
434) -> SMatrix<D, N_GROUPS, N> {
435    let mut ms = *groups;
436    ms.column_iter_mut()
437        .for_each(|mut s| s.component_mul_assign(x));
438    ms
439}
440
441fn pair_integral<D: DualNum<f64> + Copy>(mij1: D, mij2: D, eta: D, eps_ij_t: D) -> D {
442    let mut eta_i = D::one();
443    let mut j = D::zero();
444    for (ad, bd) in AD.into_iter().zip(BD) {
445        j += (eps_ij_t * (mij2 * bd[2] + mij1 * bd[1] + bd[0])
446            + (mij2 * ad[2] + mij1 * ad[1] + ad[0]))
447            * eta_i;
448        eta_i *= eta;
449    }
450    j
451}
452
453fn triplet_integral<D: DualNum<f64> + Copy>(mij1: D, mij2: D, eta: D) -> D {
454    let mut eta_i = D::one();
455    let mut j = D::zero();
456    for cd in CD {
457        j += (mij2 * cd[2] + mij1 * cd[1] + cd[0]) * eta_i;
458        eta_i *= eta;
459    }
460    j
461}
462
463#[cfg(test)]
464pub mod test {
465    use super::{
466        GcPcSaft as GcPcSaftAD, GcPcSaftParameters, ResidualHelmholtzEnergy, EPSILON_K, GROUPS, M,
467        MU, SIGMA,
468    };
469    use approx::assert_relative_eq;
470    use feos::gc_pcsaft::{GcPcSaft, GcPcSaftEosParameters, GcPcSaftRecord};
471    use feos_core::parameter::{ChemicalRecord, ParameterHetero, SegmentRecord};
472    use feos_core::{Contributions::Total, EosResult, ReferenceSystem, State};
473    use nalgebra::SVector;
474    use ndarray::arr1;
475    use quantity::{KELVIN, KILO, METER, MOL};
476    use std::collections::HashMap;
477    use std::sync::Arc;
478
479    pub fn gcpcsaft() -> EosResult<(GcPcSaftParameters<f64, 1>, Arc<GcPcSaft>)> {
480        let cr = ChemicalRecord::new(
481            Default::default(),
482            vec!["CH3".into(), ">C=O".into(), "CH2".into(), "CH3".into()],
483            None,
484        );
485        let segment_records: Vec<_> = M
486            .into_iter()
487            .zip(SIGMA)
488            .zip(EPSILON_K)
489            .zip(MU)
490            .zip(GROUPS)
491            .map(|((((m, sigma), epsilon_k), mu), g)| {
492                SegmentRecord::new(
493                    g.into(),
494                    0.0,
495                    GcPcSaftRecord::new(
496                        m,
497                        sigma,
498                        epsilon_k,
499                        Some(mu),
500                        None,
501                        None,
502                        None,
503                        None,
504                        None,
505                        None,
506                    ),
507                )
508            })
509            .collect();
510        let params = GcPcSaftEosParameters::from_segments(vec![cr], segment_records, None)?;
511        let eos = Arc::new(GcPcSaft::new(Arc::new(params)));
512        let mut groups = HashMap::new();
513        groups.insert("CH3", 2.0);
514        groups.insert(">C=O", 1.0);
515        groups.insert("CH2", 1.0);
516        let mut bonds = HashMap::new();
517        bonds.insert(["CH3", ">C=O"], 1.0);
518        bonds.insert([">C=O", "CH2"], 1.0);
519        bonds.insert(["CH2", "CH3"], 1.0);
520        let params = GcPcSaftParameters::from_groups([&groups], [&bonds]);
521        Ok((params, eos))
522    }
523
524    #[test]
525    fn test_gcpcsaft() -> EosResult<()> {
526        let (params, eos) = gcpcsaft()?;
527
528        let temperature = 300.0 * KELVIN;
529        let volume = 2.3 * METER * METER * METER;
530        let moles = arr1(&[1.3]) * KILO * MOL;
531
532        let state = State::new_nvt(&eos, temperature, volume, &moles)?;
533        let a_feos = state.residual_molar_helmholtz_energy();
534        let mu_feos = state.residual_chemical_potential();
535        let p_feos = state.pressure(Total);
536        let s_feos = state.residual_molar_entropy();
537        let h_feos = state.residual_molar_enthalpy();
538
539        let total_moles = moles.sum();
540        let t = temperature.to_reduced();
541        let v = (volume / total_moles).to_reduced();
542        let x = SVector::from_fn(|i, _| moles.get(i).convert_into(total_moles));
543        let a_ad = GcPcSaftAD::residual_molar_helmholtz_energy(&params, t, v, &x);
544        let mu_ad = GcPcSaftAD::residual_chemical_potential(&params, t, v, &x);
545        let p_ad = GcPcSaftAD::pressure(&params, t, v, &x);
546        let s_ad = GcPcSaftAD::residual_molar_entropy(&params, t, v, &x);
547        let h_ad = GcPcSaftAD::residual_molar_enthalpy(&params, t, v, &x);
548
549        println!("\nHelmholtz energy density:\n{}", a_feos.to_reduced(),);
550        println!("{a_ad}");
551        assert_relative_eq!(a_feos.to_reduced(), a_ad, max_relative = 1e-14);
552
553        println!("\nChemical potential:\n{}", mu_feos.to_reduced());
554        println!("{mu_ad}");
555        assert_relative_eq!(mu_feos.get(0).to_reduced(), mu_ad[0], max_relative = 1e-14);
556
557        println!("\nPressure:\n{}", p_feos.to_reduced());
558        println!("{p_ad}");
559        assert_relative_eq!(p_feos.to_reduced(), p_ad, max_relative = 1e-14);
560
561        println!("\nMolar entropy:\n{}", s_feos.to_reduced());
562        println!("{s_ad}");
563        assert_relative_eq!(s_feos.to_reduced(), s_ad, max_relative = 1e-14);
564
565        println!("\nMolar enthalpy:\n{}", h_feos.to_reduced());
566        println!("{h_ad}");
567        assert_relative_eq!(h_feos.to_reduced(), h_ad, max_relative = 1e-14);
568
569        Ok(())
570    }
571}
feos_ad/eos/gc_pcsaft.rs

feos_ad/eos/
gc_pcsaft.rs
