Skip to content

Add uv b3 #98

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
g-bauer merged 14 commits into from
Jan 19, 2023
Merged

Add uv b3 #98

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
21e6789
implementation uvb3-theory
Dec 22, 2022
cc4e6f6
update python.rs
Dec 23, 2022
84b2fc4
add expample, rounded parameters
Dec 23, 2022
44ab0aa
change constant names, delete unnecessary variables and comments
Jan 9, 2023
0acc89d
reduce match arms for wca
Jan 9, 2023
27e5bcd
remove multiple definition of eca one fluid properties
Jan 9, 2023
d198194
removed comments with equations
Jan 9, 2023
6a30744
add documentation for virial order
Jan 11, 2023
6f34b7f
Added generic error type to EoSError, added constructor for simple, p…
Jan 11, 2023
acc1101
cleanup error, inlcude uv-b3-theory in example, add nist simulation d…
Jan 11, 2023
590ec88
changed pub to pub(super), removed lazy statics from BH implementatio…
Jan 12, 2023
819618a
Combine example notebooks, clean up tests in mod
Jan 12, 2023
98c2c17
cleanup mod
Jan 12, 2023
e2e2889
update changelogs
g-bauer Jan 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
changed pub to pub(super), removed lazy statics from BH implementatio… 
…n, use new_simple parameters in tests
  • Loading branch information
@g-bauer
Anja Reimer authored and g-bauer committed Jan 19, 2023
commit 590ec88f54f106956cf52559927addd165530cc1
17 changes: 7 additions & 10 deletions src/uvtheory/eos/attractive_perturbation_uvb3.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,14 +151,12 @@ fn residual_virial_coefficient<D: DualNum<f64>>(p: &UVParameters, x: &Array1<D>,
let xi = x[i];

for j in 0..p.ncomponents {

let t_ij = t / p.eps_k_ij[[i, j]];
let rep_ij = p.rep_ij[[i, j]];
let att_ij = p.att_ij[[i, j]];

let q_ij = dimensionless_diameter_q_wca(t_ij, D::from(rep_ij), D::from(att_ij));


delta_b2bar +=
xi * x[j] * p.sigma_ij[[i, j]].powi(3) * delta_b2(t_ij, rep_ij, att_ij, q_ij);
}
Expand All @@ -182,10 +180,10 @@ fn residual_third_virial_coefficient<D: DualNum<f64>>(
let att_ij = p.att_ij[[i, j]];
let q_ij = dimensionless_diameter_q_wca(t_ij, D::from(rep_ij), D::from(att_ij));

// No mixing rule defined for B3 yet! The implemented rule is just taken from B2 and not correct!
let rm_ij = (rep_ij / att_ij).powd((rep_ij - att_ij).recip());
let d_ij = (d[i] / p.sigma[i] + d[j] / p.sigma[j]) * 0.5;
// Recheck mixing rule!
// No mixing rule defined for B3 yet! The implemented rule is just taken from B2 and not correct!
let rm_ij = (rep_ij / att_ij).powd((rep_ij - att_ij).recip());
let d_ij = (d[i] / p.sigma[i] + d[j] / p.sigma[j]) * 0.5;
// Recheck mixing rule!
delta_b3bar += xi
* x[j]
* p.sigma_ij[[i, j]].powi(6)
Expand Down Expand Up @@ -224,7 +222,7 @@ fn u_fraction_wca<D: DualNum<f64>>(rep_x: D, reduced_density: D, t_x: D) -> D {
+ 1.0
}

// Coefficients for IWCA
// Coefficients for IWCA
fn coefficients_wca<D: DualNum<f64>>(rep: D, att: D, d: D) -> [D; 6] {
let rep_inv = rep.recip();
let rs_x = (rep / att).powd((rep - att).recip());
Expand Down Expand Up @@ -259,14 +257,14 @@ fn coefficients_wca<D: DualNum<f64>>(rep: D, att: D, d: D) -> [D; 6] {

// Residual second virial coefficient from Revised series approximation RSAP

pub fn factorial(num: u64) -> u64 {
fn factorial(num: u64) -> u64 {
(1..=num).product()
}

fn delta_b2<D: DualNum<f64>>(reduced_temperature: D, rep: f64, att: f64, q: D) -> D {
let rm = (rep / att).powd((rep - att).recip());
let beta = reduced_temperature.recip();
let b20 = q.powi(3) * 2.0 / 3.0 * PI;
let b20 = q.powi(3) * 2.0 / 3.0 * PI;
let y = beta.exp() - 1.0;

let c1 = rep.recip() * C_B2_RSAP[0][1]
Expand Down Expand Up @@ -429,7 +427,6 @@ mod test {
let db3 = delta_b3(t_x, rm_x, rep_x, att_x, d_x, q_vdw);
assert_relative_eq!(db3.re(), -0.6591980196661884, epsilon = 1e-10);


// Full attractive perturbation:
let a = pt.helmholtz_energy(&state) / moles[0];

Expand Down
9 changes: 2 additions & 7 deletions src/uvtheory/eos/attractive_perturbation_wca.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ fn u_fraction_wca<D: DualNum<f64>>(rep_x: D, reduced_density: D) -> D {
.tanh()
}

pub fn one_fluid_properties<D: DualNum<f64>>(
pub(super) fn one_fluid_properties<D: DualNum<f64>>(
p: &UVParameters,
x: &Array1<D>,
t: D,
Expand Down Expand Up @@ -315,16 +315,13 @@ mod test {
let q_vdw = dimensionless_diameter_q_wca(t_x, rep_x, att_x);
let b21u = delta_b12u(t_x, mean_field_constant_x, weighted_sigma3_ij, q_vdw, rm_x)
/ p.sigma[0].powi(3);
//assert!(b21u.re() == -1.02233216);
assert_relative_eq!(b21u.re(), -1.02233215790525, epsilon = 1e-12);

let i_wca =
correlation_integral_wca(rho_x, mean_field_constant_x, rep_x, att_x, d_x, q_vdw, rm_x);

let delta_a1u = state.partial_density.sum() / t_x * i_wca * 2.0 * PI * weighted_sigma3_ij;

//dbg!(delta_a1u);
//assert!(delta_a1u.re() == -1.1470186919354);
assert_relative_eq!(delta_a1u.re(), -1.52406840346272, epsilon = 1e-6);

let u_fraction_wca = u_fraction_wca(
Expand All @@ -336,10 +333,8 @@ mod test {
dbg!(b2bar);
assert_relative_eq!(b2bar.re(), -1.09102560732964, epsilon = 1e-12);
dbg!(u_fraction_wca);
//assert!(u_fraction_WCA.re() == 0.743451055308332);
assert_relative_eq!(u_fraction_wca.re(), 0.997069754340431, epsilon = 1e-5);

//assert!(b2bar.re() == -1.00533412744652);
assert_relative_eq!(u_fraction_wca.re(), 0.997069754340431, epsilon = 1e-5);

let a_test = delta_a1u
+ (-u_fraction_wca + 1.0)
Expand Down
91 changes: 36 additions & 55 deletions src/uvtheory/eos/hard_sphere_bh.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,24 +1,22 @@
use crate::uvtheory::parameters::UVParameters;
use feos_core::{HelmholtzEnergyDual, StateHD};
use lazy_static::lazy_static;
use ndarray::prelude::*;
use num_dual::DualNum;
use std::fmt;
use std::sync::Arc;

lazy_static! {
static ref BH_CONSTANTS_ETA_B: Array2<f64> = arr2(&[
[-0.960919783, -0.921097447],
[-0.547468020, -3.508014069],
[-2.253750186, 3.581161364],
]);
static ref BH_CONSTANTS_ETA_A: Array2<f64> = arr2(&[
[-1.217417282, 6.754987582, -0.5919326153, -28.99719604],
[1.579548775, -26.93879416, 0.3998915410, 106.9446266],
[-1.993990512, 44.11863355, -40.10916106, -29.6130848],
[0.0, 0.0, 0.0, 0.0],
]);
}
const BH_CONSTANTS_ETA_B: [[f64; 2]; 3] = [
[-0.960919783, -0.921097447],
[-0.547468020, -3.508014069],
[-2.253750186, 3.581161364],
];

const BH_CONSTANTS_ETA_A: [[f64; 4]; 4] = [
[-1.217417282, 6.754987582, -0.5919326153, -28.99719604],
[1.579548775, -26.93879416, 0.3998915410, 106.9446266],
[-1.993990512, 44.11863355, -40.10916106, -29.6130848],
[0.0, 0.0, 0.0, 0.0],
];

#[derive(Debug, Clone)]
pub struct HardSphereBH {
Expand Down Expand Up @@ -48,7 +46,7 @@ impl fmt::Display for HardSphereBH {
/// Dimensionless Hard-sphere diameter according to Barker-Henderson division.
/// Eq. S23 and S24.
///
pub fn diameter_bh<D: DualNum<f64>>(parameters: &UVParameters, temperature: D) -> Array1<D> {
pub(super) fn diameter_bh<D: DualNum<f64>>(parameters: &UVParameters, temperature: D) -> Array1<D> {
parameters
.cd_bh_pure
.iter()
Expand All @@ -63,7 +61,7 @@ pub fn diameter_bh<D: DualNum<f64>>(parameters: &UVParameters, temperature: D) -
.collect()
}

pub fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> [D; 4] {
pub(super) fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> [D; 4] {
let mut zeta: [D; 4] = [D::zero(), D::zero(), D::zero(), D::zero()];
for i in 0..partial_density.len() {
for k in 0..4 {
Expand All @@ -74,13 +72,16 @@ pub fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>)
zeta
}

pub fn packing_fraction<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> D {
pub(super) fn packing_fraction<D: DualNum<f64>>(
partial_density: &Array1<D>,
diameter: &Array1<D>,
) -> D {
(0..partial_density.len()).fold(D::zero(), |acc, i| {
acc + partial_density[i] * diameter[i].powi(3) * (std::f64::consts::PI / 6.0)
})
}

pub fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) -> D {
pub(super) fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) -> D {
let mut zeta: [D; 2] = [D::zero(), D::zero()];
for i in 0..molefracs.len() {
for k in 0..2 {
Expand All @@ -90,7 +91,7 @@ pub fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) ->
zeta[0] / zeta[1]
}

pub fn packing_fraction_b<D: DualNum<f64>>(
pub(super) fn packing_fraction_b<D: DualNum<f64>>(
parameters: &UVParameters,
diameter: &Array1<D>,
eta: D,
Expand All @@ -100,36 +101,39 @@ pub fn packing_fraction_b<D: DualNum<f64>>(
let tau =
-(diameter[i] / parameters.sigma[i] + diameter[j] / parameters.sigma[j]) * 0.5 + 1.0; //dimensionless
let tau2 = tau * tau;

let c = arr1(&[
tau * BH_CONSTANTS_ETA_B[[0, 0]] + tau2 * BH_CONSTANTS_ETA_B[[0, 1]],
tau * BH_CONSTANTS_ETA_B[[1, 0]] + tau2 * BH_CONSTANTS_ETA_B[[1, 1]],
tau * BH_CONSTANTS_ETA_B[[2, 0]] + tau2 * BH_CONSTANTS_ETA_B[[2, 1]],
tau * BH_CONSTANTS_ETA_B[0][0] + tau2 * BH_CONSTANTS_ETA_B[0][1],
tau * BH_CONSTANTS_ETA_B[1][0] + tau2 * BH_CONSTANTS_ETA_B[1][1],
tau * BH_CONSTANTS_ETA_B[2][0] + tau2 * BH_CONSTANTS_ETA_B[2][1],
]);
eta + eta * c[0] + eta * eta * c[1] + eta.powi(3) * c[2]
})
}

pub fn packing_fraction_a<D: DualNum<f64>>(
pub(super) fn packing_fraction_a<D: DualNum<f64>>(
parameters: &UVParameters,
diameter: &Array1<D>,
eta: D,
) -> Array2<D> {
let n = parameters.att.len();
Array2::from_shape_fn((n, n), |(i, j)| {
let tau =
-(diameter[i] / parameters.sigma[i] + diameter[j] / parameters.sigma[j]) * 0.5 + 1.0; //dimensionless//-(diameter[i] + diameter[j]) * 0.5 + 1.0;
-(diameter[i] / parameters.sigma[i] + diameter[j] / parameters.sigma[j]) * 0.5 + 1.0;
let tau2 = tau * tau;
let rep_inv = 1.0 / parameters.rep_ij[[i, j]];

let c = arr1(&[
tau * (BH_CONSTANTS_ETA_A[[0, 0]] + BH_CONSTANTS_ETA_A[[0, 1]] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[[0, 2]] + BH_CONSTANTS_ETA_A[[0, 3]] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[[1, 0]] + BH_CONSTANTS_ETA_A[[1, 1]] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[[1, 2]] + BH_CONSTANTS_ETA_A[[1, 3]] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[[2, 0]] + BH_CONSTANTS_ETA_A[[2, 1]] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[[2, 2]] + BH_CONSTANTS_ETA_A[[2, 3]] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[[3, 0]] + BH_CONSTANTS_ETA_A[[3, 1]] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[[3, 2]] + BH_CONSTANTS_ETA_A[[3, 3]] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[0][0] + BH_CONSTANTS_ETA_A[0][1] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[0][2] + BH_CONSTANTS_ETA_A[0][3] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[1][0] + BH_CONSTANTS_ETA_A[1][1] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[1][2] + BH_CONSTANTS_ETA_A[1][3] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[2][0] + BH_CONSTANTS_ETA_A[2][1] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[2][2] + BH_CONSTANTS_ETA_A[2][3] * rep_inv),
tau * (BH_CONSTANTS_ETA_A[3][0] + BH_CONSTANTS_ETA_A[3][1] * rep_inv)
+ tau2 * (BH_CONSTANTS_ETA_A[3][2] + BH_CONSTANTS_ETA_A[3][3] * rep_inv),
]);

eta + eta * c[0] + eta * eta * c[1] + eta.powi(3) * c[2] + eta.powi(4) * c[3]
})
}
Expand Down Expand Up @@ -158,27 +162,4 @@ mod test {
0.95583586434435486
);
}

// #[test]
// fn helmholtz_energy() {
// let p = methane_parameters(12);
// let p = test_parameters(12, 1.0, 1.0);
// let reduced_density = 1.0;
// let reduced_temperature = 4.0;

// let hs = HardSphere {
// parameters: Arc::new(p.clone()),
// };
// let particles = arr1(&[1000.0]);
// let n = &particles / 6.02214076e23;
// let volume = particles[0] / reduced_density * p.sigma[0].powi(3);
// let temperature = reduced_temperature * p.epsilon_k[0];
// let s = StateHD::new(temperature, volume, n.clone());
// dbg!(particles[0] / volume * p.sigma[0].powi(3));
// dbg!(&s.temperature);
// dbg!(&s.volume);
// dbg!(&s.moles);
// let a = hs.helmholtz_energy(&s);
// assert_relative_eq!(a / particles[0], 1.9859860794723039, epsilon = 1e-10)
// }
}
32 changes: 19 additions & 13 deletions src/uvtheory/eos/hard_sphere_wca.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ use std::f64::consts::PI;
use std::fmt;
use std::sync::Arc;

pub const WCA_CONSTANTS_Q: [[f64; 4]; 3] = [
pub(super) const WCA_CONSTANTS_Q: [[f64; 4]; 3] = [
[1.92840364363978, 4.43165896265079E-01, 0.0, 0.0],
[
5.20120816141761E-01,
Expand Down Expand Up @@ -36,14 +36,14 @@ const WCA_CONSTANTS_ETA_B: [[f64; 2]; 3] = [
];

// New Fit to numerical integrals for uv-B3-theory
pub const WCA_CONSTANTS_ETA_A_UVB3: [[f64; 4]; 4] = [
pub(super) const WCA_CONSTANTS_ETA_A_UVB3: [[f64; 4]; 4] = [
[2.64043218, -1.2184421, -22.90786387, 0.96433414],
[-16.75643936, 30.83929771, 73.08711814, -166.57701616],
[19.53170162, -88.87955657, -76.51387192, 443.68942745],
[-3.77740877, 83.04694547, 21.62502721, -304.8643176],
];

pub const WCA_CONSTANTS_ETA_B_UVB3: [[f64; 2]; 3] = [
pub(super) const WCA_CONSTANTS_ETA_B_UVB3: [[f64; 2]; 3] = [
[2.19821588, -20.45005484],
[-13.47050687, 56.65701375],
[12.90119266, -42.71680606],
Expand Down Expand Up @@ -74,7 +74,10 @@ impl fmt::Display for HardSphereWCA {
}

/// Dimensionless Hard-sphere diameter according to Weeks-Chandler-Andersen division.
pub fn diameter_wca<D: DualNum<f64>>(parameters: &UVParameters, temperature: D) -> Array1<D> {
pub(super) fn diameter_wca<D: DualNum<f64>>(
parameters: &UVParameters,
temperature: D,
) -> Array1<D> {
parameters
.sigma
.iter()
Expand All @@ -93,7 +96,7 @@ pub fn diameter_wca<D: DualNum<f64>>(parameters: &UVParameters, temperature: D)
.collect()
}

pub fn dimensionless_diameter_q_wca<D: DualNum<f64>>(t_x: D, rep_x: D, att_x: D) -> D {
pub(super) fn dimensionless_diameter_q_wca<D: DualNum<f64>>(t_x: D, rep_x: D, att_x: D) -> D {
let nu = rep_x;
let n = att_x;
let rs = (nu / n).powd((nu - n).recip());
Expand All @@ -119,7 +122,7 @@ pub fn dimensionless_diameter_q_wca<D: DualNum<f64>>(t_x: D, rep_x: D, att_x: D)
* rs
}

pub fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> [D; 4] {
pub(super) fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> [D; 4] {
let mut zeta: [D; 4] = [D::zero(), D::zero(), D::zero(), D::zero()];
for i in 0..partial_density.len() {
for k in 0..4 {
Expand All @@ -130,13 +133,16 @@ pub fn zeta<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>)
zeta
}

pub fn packing_fraction<D: DualNum<f64>>(partial_density: &Array1<D>, diameter: &Array1<D>) -> D {
pub(super) fn packing_fraction<D: DualNum<f64>>(
partial_density: &Array1<D>,
diameter: &Array1<D>,
) -> D {
(0..partial_density.len()).fold(D::zero(), |acc, i| {
acc + partial_density[i] * diameter[i].powi(3) * (std::f64::consts::PI / 6.0)
})
}

pub fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) -> D {
pub(super) fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) -> D {
let mut zeta: [D; 2] = [D::zero(), D::zero()];
for i in 0..molefracs.len() {
for k in 0..2 {
Expand All @@ -147,7 +153,7 @@ pub fn zeta_23<D: DualNum<f64>>(molefracs: &Array1<D>, diameter: &Array1<D>) ->
}

#[inline]
pub fn dimensionless_length_scale<D: DualNum<f64>>(
pub(super) fn dimensionless_length_scale<D: DualNum<f64>>(
parameters: &UVParameters,
temperature: D,
) -> Array1<D> {
Expand All @@ -166,7 +172,7 @@ pub fn dimensionless_length_scale<D: DualNum<f64>>(

#[inline]

pub fn packing_fraction_b<D: DualNum<f64>>(
pub(super) fn packing_fraction_b<D: DualNum<f64>>(
parameters: &UVParameters,
eta: D,
temperature: D,
Expand All @@ -188,7 +194,7 @@ pub fn packing_fraction_b<D: DualNum<f64>>(
})
}

pub fn packing_fraction_b_uvb3<D: DualNum<f64>>(
pub(super) fn packing_fraction_b_uvb3<D: DualNum<f64>>(
parameters: &UVParameters,
eta: D,
temperature: D,
Expand All @@ -210,7 +216,7 @@ pub fn packing_fraction_b_uvb3<D: DualNum<f64>>(
})
}

pub fn packing_fraction_a<D: DualNum<f64>>(
pub(super) fn packing_fraction_a<D: DualNum<f64>>(
parameters: &UVParameters,
eta: D,
temperature: D,
Expand Down Expand Up @@ -239,7 +245,7 @@ pub fn packing_fraction_a<D: DualNum<f64>>(
})
}

pub fn packing_fraction_a_uvb3<D: DualNum<f64>>(
pub(super) fn packing_fraction_a_uvb3<D: DualNum<f64>>(
parameters: &UVParameters,
eta: D,
temperature: D,
Expand Down