#![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>> {

// correct for different k_B in old code

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>> {

// correct for different k_B in old code

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>> {

// correct for different k_B in old code

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

}