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

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
removed comments with equations
  • Loading branch information
@g-bauer
Anja Reimer authored and g-bauer committed Jan 19, 2023
commit d1981948df0c3e7c0958422482215fac8fcb9847
28 changes: 11 additions & 17 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 @@ -91,7 +91,7 @@ impl fmt::Display for AttractivePerturbationUVB3 {
}

impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for AttractivePerturbationUVB3 {
/// Helmholtz energy for attractive perturbation, eq. 52
/// Helmholtz energy for attractive perturbation
fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
let p = &self.parameters;
let t = state.temperature;
Expand All @@ -111,7 +111,6 @@ impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for AttractivePerturbationUVB3 {

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

// state.partial_density.sum() / t_x * i_wca * 2.0 * PI * weighted_sigma3_ij;
let u_fraction_wca = u_fraction_wca(
rep_x,
density * (x * &p.sigma.mapv(|s| s.powi(3))).sum(),
Expand All @@ -132,7 +131,6 @@ impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for AttractivePerturbationUVB3 {
}
}

// (S43) & (S53)
fn delta_b12u<D: DualNum<f64>>(
t_x: D,
mean_field_constant_x: D,
Expand All @@ -153,14 +151,14 @@ fn residual_virial_coefficient<D: DualNum<f64>>(p: &UVParameters, x: &Array1<D>,
let xi = x[i];

for j in 0..p.ncomponents {
//let q_ij = (q[i] / p.sigma[i] + q[j] / p.sigma[j]) * 0.5;

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));

// Recheck mixing rule!

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

for j in 0..p.ncomponents {
//let q_ij = (q[i] / p.sigma[i] + q[j] / p.sigma[j]) * 0.5;
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));

// Recheck mixing rule!
let rm_ij = (rep_ij / att_ij).powd((rep_ij - att_ij).recip()); // Anja frei erfunden
let d_ij = (d[i] / p.sigma[i] + d[j] / p.sigma[j]) * 0.5; // Recheck mixing rule!
// // 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 All @@ -215,7 +212,6 @@ fn correlation_integral_wca<D: DualNum<f64>>(

/// U-fraction with low temperature correction omega
fn u_fraction_wca<D: DualNum<f64>>(rep_x: D, reduced_density: D, t_x: D) -> D {
// let omega = (-t_x.powi(2) * rep_x * CU_WCA[5]).exp() * (t_x.recip() * CU_WCA[6]).exp();
let omega = if t_x.re() < 175.0 {
(-t_x * CU_WCA[5] * (reduced_density - CU_WCA[6]).powi(2)).exp()
* ((t_x * CU_WCA[7]).tanh().recip() - 1.0).powi(2)
Expand All @@ -228,7 +224,7 @@ fn u_fraction_wca<D: DualNum<f64>>(rep_x: D, reduced_density: D, t_x: D) -> D {
+ 1.0
}

// Coefficients for IWCA from eq. (S55)
// 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 @@ -268,9 +264,9 @@ pub fn factorial(num: u64) -> u64 {
}

fn delta_b2<D: DualNum<f64>>(reduced_temperature: D, rep: f64, att: f64, q: D) -> D {
let rm = (rep / att).powd((rep - att).recip()); // Check mixing rule!!
let rm = (rep / att).powd((rep - att).recip());
let beta = reduced_temperature.recip();
let b20 = q.powi(3) * 2.0 / 3.0 * PI; // eq. (16)
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 @@ -353,7 +349,7 @@ fn delta_b3<D: DualNum<f64>>(t_x: D, rm_x: f64, rep_x: f64, _att_x: f64, d_x: D,

let b3 = b30 + b31 * beta + b32 * beta.powi(2) + b33 * beta.powi(3) + b34 * beta.powi(4);

// Watch out: Not implemented for mixtures!!!!
// Watch out: Not defined for mixtures!
let tau = -d_x + rm_x;
let tau2 = tau * tau;
let rep_inv = rep_x.recip();
Expand Down Expand Up @@ -433,8 +429,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);

// let db31 =
// delta_b31u(t_x, weighted_sigma3_ij, rm_x, rep_x, att_x, d_x) / p.sigma[0].powi(6);

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