#![allow(clippy::excessive_precision)]
use approx::assert_relative_eq;
use feos_core::parameter::{IdentifierOption, Parameter};
use feos_core::{Contributions, PhaseEquilibrium, State};
use feos_dft::fundamental_measure_theory::FMTVersion;
use feos_dft::interface::PlanarInterface;
use feos_pcsaft::{PcSaft, PcSaftFunctional, PcSaftParameters};
use ndarray::{arr1, Axis};
use quantity::si::*;
use std::error::Error;
use std::rc::Rc;
#[test]
#[allow(non_snake_case)]
fn test_bulk_implementations() -> Result<(), Box<dyn Error>> {
let KB_old = 1.38064852e-23 * JOULE / KELVIN;
let NAV_old = 6.022140857e23 / MOL;
let params = Rc::new(PcSaftParameters::from_json(
vec!["water_np"],
"tests/test_parameters.json",
None,
IdentifierOption::Name,
)?);
let eos = Rc::new(PcSaft::new(params.clone()));
let func_pure = Rc::new(PcSaftFunctional::new(params.clone()));
let func_full = Rc::new(PcSaftFunctional::new_full(
params,
FMTVersion::KierlikRosinberg,
));
let t = 300.0 * KELVIN;
let v = 0.002 * METER.powi(3) * NAV / NAV_old;
let n = arr1(&[1.5]) * MOL;
let state = State::new_nvt(&eos, t, v, &n)?;
let state_pure = State::new_nvt(&func_pure, t, v, &n)?;
let state_full = State::new_nvt(&func_full, t, v, &n)?;
let p = state.pressure_contributions();
let p_pure = state_pure.pressure_contributions();
let p_full = state_full.pressure_contributions();
println!("{}: {}", p[0].0, p[0].1);
println!("{}: {}", p_pure[0].0, p_pure[0].1);
println!("{}: {}", p_full[0].0, p_full[0].1);
println!();
println!("{:20}: {}", p[1].0, p[1].1);
println!("{:20}: {}", p_pure[1].0, p_pure[1].1);
println!("{:20}: {}", p_full[1].0, p_full[1].1);
println!();
println!("{:21}: {}", p[2].0, p[2].1);
println!("{:21}: {}", p_pure[2].0, p_pure[2].1);
println!("{:21}: {}", p_full[2].0, p_full[2].1);
println!();
println!("{:21}: {}", p[3].0, p[3].1);
println!("{:21}: {}", p_pure[3].0, p_pure[3].1 + p_pure[4].1);
println!("{:21}: {}", p_full[3].0, p_full[3].1 + p_full[5].1);
println!();
println!("{:21}: {}", p[4].0, p[4].1);
println!("{:21}: {}", p_full[4].0, p_full[4].1);
let ideal_gas = 1.8707534688259309 * MEGA * PASCAL * KB / KB_old;
assert_relative_eq!(p[0].1, ideal_gas, max_relative = 1e-14,);
assert_relative_eq!(p_pure[0].1, ideal_gas, max_relative = 1e-14,);
assert_relative_eq!(p_full[0].1, ideal_gas, max_relative = 1e-14,);
let hard_sphere = 54.7102253882827583 * KILO * PASCAL * KB / KB_old;
assert_relative_eq!(p[1].1, hard_sphere, max_relative = 1e-14,);
assert_relative_eq!(p_full[1].1, hard_sphere, max_relative = 1e-14,);
let hard_chains = -2.0847750028499230 * KILO * PASCAL * KB / KB_old;
assert_relative_eq!(p[2].1, hard_chains, max_relative = 1e-14,);
assert_relative_eq!(p_pure[2].1 + p_pure[4].1, hard_chains, max_relative = 2e-13,);
assert_relative_eq!(p_full[2].1 + p_full[5].1, hard_chains, max_relative = 2e-13,);
let dispersion = -262.895932352779993 * KILO * PASCAL * KB / KB_old;
assert_relative_eq!(p[3].1, dispersion, max_relative = 1e-14,);
assert_relative_eq!(p_pure[3].1, dispersion, max_relative = 1e-14,);
assert_relative_eq!(p_full[3].1, dispersion, max_relative = 1e-14,);
let association = -918.3899928262694630 * KILO * PASCAL * KB / KB_old;
assert_relative_eq!(p[4].1, association, max_relative = 1e-14,);
assert_relative_eq!(p_full[4].1, association, max_relative = 1e-14,);
assert_relative_eq!(p_pure[1].1, hard_sphere + association, max_relative = 1e-14,);
Ok(())
}
#[test]
#[allow(non_snake_case)]
fn test_dft_propane() -> Result<(), Box<dyn Error>> {
let KB_old = 1.38064852e-23 * JOULE / KELVIN;
let NAV_old = 6.022140857e23 / MOL;
let params = Rc::new(PcSaftParameters::from_json(
vec!["propane"],
"tests/test_parameters.json",
None,
IdentifierOption::Name,
)?);
let func_pure = Rc::new(PcSaftFunctional::new(params.clone()));
let func_full = Rc::new(PcSaftFunctional::new_full(
params.clone(),
FMTVersion::KierlikRosinberg,
));
let func_full_vec = Rc::new(PcSaftFunctional::new_full(params, FMTVersion::WhiteBear));
let t = 200.0 * KELVIN;
let w = 150.0 * ANGSTROM;
let points = 2048;
let tc = State::critical_point(&func_pure, None, None, Default::default())?.temperature;
let vle_pure = PhaseEquilibrium::pure(&func_pure, t, None, Default::default())?;
let vle_full = PhaseEquilibrium::pure(&func_full, t, None, Default::default())?;
let vle_full_vec = PhaseEquilibrium::pure(&func_full_vec, t, None, Default::default())?;
let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc)?.solve(None)?;
let profile_full = PlanarInterface::from_tanh(&vle_full, points, w, tc)?.solve(None)?;
let profile_full_vec = PlanarInterface::from_tanh(&vle_full_vec, points, w, tc)?.solve(None)?;
let _ = func_pure.solve_pdgt(&vle_pure, 198, 0, None)?;
println!(
"pure {} {} {} {}",
profile_pure.surface_tension.unwrap(),
vle_pure.vapor().density,
vle_pure.liquid().density,
func_pure.solve_pdgt(&vle_pure, 198, 0, None)?.1
);
println!(
"full {} {} {} {}",
profile_full.surface_tension.unwrap(),
vle_full.vapor().density,
vle_full.liquid().density,
func_full.solve_pdgt(&vle_full, 198, 0, None)?.1
);
println!(
"vec {} {} {} {}",
profile_full_vec.surface_tension.unwrap(),
vle_full_vec.vapor().density,
vle_full_vec.liquid().density,
func_full_vec.solve_pdgt(&vle_full_vec, 198, 0, None)?.1
);
let vapor_density = 12.2557486248527745 * MOL / METER.powi(3) * NAV_old / NAV;
assert_relative_eq!(
vle_pure.vapor().density,
vapor_density,
max_relative = 1e-13,
);
assert_relative_eq!(
vle_full.vapor().density,
vapor_density,
max_relative = 1e-13,
);
assert_relative_eq!(
vle_full_vec.vapor().density,
vapor_density,
max_relative = 1e-13,
);
let liquid_density = 13.8941749145544549 * KILO * MOL / METER.powi(3) * NAV_old / NAV;
assert_relative_eq!(
vle_pure.liquid().density,
liquid_density,
max_relative = 1e-13,
);
assert_relative_eq!(
vle_full.liquid().density,
liquid_density,
max_relative = 1e-13,
);
assert_relative_eq!(
vle_full_vec.liquid().density,
liquid_density,
max_relative = 1e-13,
);
let surface_tension = 19.9931025166113692 * MILLI * NEWTON / METER * KB / KB_old;
let surface_tension_kr = 19.9863313312996169 * MILLI * NEWTON / METER * KB / KB_old;
assert_relative_eq!(
profile_pure.surface_tension.unwrap(),
surface_tension,
max_relative = 1e-4,
);
assert_relative_eq!(
profile_full.surface_tension.unwrap(),
surface_tension_kr,
max_relative = 1e-4,
);
assert_relative_eq!(
profile_full_vec.surface_tension.unwrap(),
surface_tension,
max_relative = 1e-4,
);
let surface_tension_pdgt = 20.2849756479219039 * MILLI * NEWTON / METER * KB / KB_old;
let surface_tension_pdgt_kr = 20.2785079953823342 * MILLI * NEWTON / METER * KB / KB_old;
assert_relative_eq!(
func_pure.solve_pdgt(&vle_pure, 198, 0, None)?.1,
surface_tension_pdgt,
max_relative = 1e-10,
);
assert_relative_eq!(
func_full.solve_pdgt(&vle_full, 198, 0, None)?.1,
surface_tension_pdgt_kr,
max_relative = 1e-10,
);
assert_relative_eq!(
func_full_vec.solve_pdgt(&vle_full_vec, 198, 0, None)?.1,
surface_tension_pdgt,
max_relative = 1e-10,
);
Ok(())
}
#[test]
#[allow(non_snake_case)]
fn test_dft_water() -> Result<(), Box<dyn Error>> {
let KB_old = 1.38064852e-23 * JOULE / KELVIN;
let NAV_old = 6.022140857e23 / MOL;
let params = Rc::new(PcSaftParameters::from_json(
vec!["water_np"],
"tests/test_parameters.json",
None,
IdentifierOption::Name,
)?);
let func_pure = Rc::new(PcSaftFunctional::new(params.clone()));
let func_full_vec = Rc::new(PcSaftFunctional::new_full(params, FMTVersion::WhiteBear));
let t = 400.0 * KELVIN;
let w = 120.0 * ANGSTROM;
let points = 2048;
let tc = State::critical_point(&func_pure, None, None, Default::default())?.temperature;
let vle_pure = PhaseEquilibrium::pure(&func_pure, t, None, Default::default())?;
let vle_full_vec = PhaseEquilibrium::pure(&func_full_vec, t, None, Default::default())?;
let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc)?.solve(None)?;
let profile_full_vec = PlanarInterface::from_tanh(&vle_full_vec, points, w, tc)?.solve(None)?;
println!(
"pure {} {} {}",
profile_pure.surface_tension.unwrap(),
vle_pure.vapor().density,
vle_pure.liquid().density
);
println!(
"vec {} {} {}",
profile_full_vec.surface_tension.unwrap(),
vle_full_vec.vapor().density,
vle_full_vec.liquid().density
);
let vapor_density = 75.8045715345905222 * MOL / METER.powi(3) * NAV_old / NAV;
assert_relative_eq!(
vle_pure.vapor().density,
vapor_density,
max_relative = 1e-13,
);
assert_relative_eq!(
vle_full_vec.vapor().density,
vapor_density,
max_relative = 1e-13,
);
let liquid_density = 47.8480850281608454 * KILO * MOL / METER.powi(3) * NAV_old / NAV;
assert_relative_eq!(
vle_pure.liquid().density,
liquid_density,
max_relative = 1e-13,
);
assert_relative_eq!(
vle_full_vec.liquid().density,
liquid_density,
max_relative = 1e-13,
);
let surface_tension = 70.1419567481980408 * MILLI * NEWTON / METER * KB / KB_old;
assert_relative_eq!(
profile_pure.surface_tension.unwrap(),
surface_tension,
max_relative = 1e-4,
);
assert_relative_eq!(
profile_full_vec.surface_tension.unwrap(),
surface_tension,
max_relative = 1e-4,
);
let surface_tension_pdgt = 72.9496195135527188 * MILLI * NEWTON / METER * KB / KB_old;
assert_relative_eq!(
func_pure.solve_pdgt(&vle_pure, 198, 0, None)?.1,
surface_tension_pdgt,
max_relative = 1e-10,
);
assert_relative_eq!(
func_full_vec.solve_pdgt(&vle_full_vec, 198, 0, None)?.1,
surface_tension_pdgt,
max_relative = 1e-10,
);
Ok(())
}
#[test]
fn test_entropy_bulk_values() -> Result<(), Box<dyn Error>> {
let params = PcSaftParameters::from_json(
vec!["water_np"],
"tests/test_parameters.json",
None,
IdentifierOption::Name,
)?;
let func = Rc::new(PcSaftFunctional::new(Rc::new(params)));
let vle = PhaseEquilibrium::pure(&func, 350.0 * KELVIN, None, Default::default())?;
let profile = PlanarInterface::from_pdgt(&vle, 2048)?.solve(None)?;
let s_res = profile
.profile
.entropy_density(Contributions::ResidualNvt)?;
let s_tot = profile.profile.entropy_density(Contributions::Total)?;
println!(
"Density:\n{}",
profile.profile.density.index_axis(Axis(0), 0).to_owned()
);
println!(
"liquid: {}, vapor: {}",
profile.vle.liquid().density,
profile.vle.vapor().density
);
println!("\nResidual:\n{}", s_res);
println!(
"liquid: {}, vapor: {}",
profile.vle.liquid().entropy(Contributions::ResidualNvt) / profile.vle.liquid().volume,
profile.vle.vapor().entropy(Contributions::ResidualNvt) / profile.vle.vapor().volume
);
println!("\nTotal:\n{}", s_tot);
println!(
"liquid: {}, vapor: {}",
profile.vle.liquid().entropy(Contributions::Total) / profile.vle.liquid().volume,
profile.vle.vapor().entropy(Contributions::Total) / profile.vle.vapor().volume
);
assert_relative_eq!(
s_res.get(0),
profile.vle.liquid().entropy(Contributions::ResidualNvt) / profile.vle.liquid().volume,
max_relative = 1e-8,
);
assert_relative_eq!(
s_res.get(2047),
profile.vle.vapor().entropy(Contributions::ResidualNvt) / profile.vle.vapor().volume,
max_relative = 1e-8,
);
assert_relative_eq!(
s_tot.get(0),
profile.vle.liquid().entropy(Contributions::Total) / profile.vle.liquid().volume,
max_relative = 1e-8,
);
assert_relative_eq!(
s_tot.get(2047),
profile.vle.vapor().entropy(Contributions::Total) / profile.vle.vapor().volume,
max_relative = 1e-8,
);
Ok(())
}