#[cfg(feature = "experimental-randnla")]
mod bench_randnla_vs_eht {
use sci_form::beta::rand_nla::*;
use sci_form::eht::solve_eht;
use std::time::Instant;
fn water() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![8, 1, 1],
vec![
[0.0, 0.0, 0.117],
[0.0, 0.757, -0.469],
[0.0, -0.757, -0.469],
],
)
}
fn ethanol() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![6, 6, 8, 1, 1, 1, 1, 1, 1],
vec![
[0.0, 0.0, 0.0],
[1.52, 0.0, 0.0],
[2.13, 1.21, 0.0],
[-0.36, 1.03, 0.0],
[-0.36, -0.51, 0.89],
[-0.36, -0.51, -0.89],
[1.88, -0.51, 0.89],
[1.88, -0.51, -0.89],
[3.09, 1.21, 0.0],
],
)
}
#[test]
fn b1_eigenvalue_deviation_water() {
let (elems, pos) = water();
let eht = solve_eht(&elems, &pos, None).unwrap();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let config = RandNlaConfig {
seed: 42,
..Default::default()
};
let (rand_energies, _, info) = solve_eht_randnla(&h, &s, &config);
let n = eht.energies.len().min(rand_energies.len());
let mut max_dev = 0.0_f64;
for i in 0..n {
let dev = (eht.energies[i] - rand_energies[i]).abs();
max_dev = max_dev.max(dev);
}
println!("[B1] Water eigenvalue max deviation: {:.6} eV", max_dev);
println!("[B1] RandNLA residual error: {:.2e}", info.residual_error);
println!("[B1] Used fallback: {}", info.used_fallback);
assert!(
max_dev < 0.1 || info.used_fallback,
"Eigenvalue deviation {:.4} eV too large without fallback",
max_dev
);
}
#[test]
fn b1_homo_lumo_gap_agreement() {
let (elems, pos) = ethanol();
let eht = solve_eht(&elems, &pos, None).unwrap();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let config = RandNlaConfig::default();
let (rand_energies, _, _) = solve_eht_randnla(&h, &s, &config);
let homo = eht.homo_index;
let lumo = eht.lumo_index;
let eht_gap = eht.gap;
let rand_gap = if lumo < rand_energies.len() && homo < rand_energies.len() {
rand_energies[lumo] - rand_energies[homo]
} else {
0.0
};
let gap_dev = (eht_gap - rand_gap).abs();
println!(
"[B1] Ethanol HOMO-LUMO gap: EHT={:.3} eV, RandNLA={:.3} eV, Δ={:.4} eV",
eht_gap, rand_gap, gap_dev
);
assert!(
gap_dev < 0.5,
"Gap deviation {:.3} eV exceeds tolerance",
gap_dev
);
}
#[test]
fn b1_timing_comparison() {
let (elems, pos) = ethanol();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let t0 = Instant::now();
for _ in 0..50 {
let _ = solve_eht(&elems, &pos, None);
}
let eht_time = t0.elapsed().as_micros() as f64 / 50.0;
let config = RandNlaConfig::default();
let t0 = Instant::now();
for _ in 0..50 {
let _ = solve_eht_randnla(&h, &s, &config);
}
let rand_time = t0.elapsed().as_micros() as f64 / 50.0;
println!(
"[B1] Timing: EHT={:.0} µs, RandNLA={:.0} µs, ratio={:.2}x",
eht_time,
rand_time,
rand_time / eht_time.max(1.0)
);
}
}
#[cfg(feature = "experimental-kpm")]
mod bench_kpm_vs_eht {
use sci_form::beta::kpm::*;
use sci_form::eht::solve_eht;
use std::time::Instant;
fn water() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![8, 1, 1],
vec![
[0.0, 0.0, 0.117],
[0.0, 0.757, -0.469],
[0.0, -0.757, -0.469],
],
)
}
fn ethanol() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![6, 6, 8, 1, 1, 1, 1, 1, 1],
vec![
[0.0, 0.0, 0.0],
[1.52, 0.0, 0.0],
[2.13, 1.21, 0.0],
[-0.36, 1.03, 0.0],
[-0.36, -0.51, 0.89],
[-0.36, -0.51, -0.89],
[1.88, -0.51, 0.89],
[1.88, -0.51, -0.89],
[3.09, 1.21, 0.0],
],
)
}
#[test]
fn b2_dos_spectral_bounds_agree() {
let (elems, pos) = water();
let eht = solve_eht(&elems, &pos, None).unwrap();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let (kpm_emin, kpm_emax) = estimate_spectral_bounds(&h);
let eht_emin = eht.energies.iter().cloned().fold(f64::INFINITY, f64::min);
let eht_emax = eht
.energies
.iter()
.cloned()
.fold(f64::NEG_INFINITY, f64::max);
println!(
"[B2] EHT spectral range: [{:.2}, {:.2}] eV",
eht_emin, eht_emax
);
println!(
"[B2] KPM spectral range: [{:.2}, {:.2}] eV",
kpm_emin, kpm_emax
);
assert!(kpm_emin <= eht_emin + 1.0, "KPM e_min too high");
assert!(kpm_emax >= eht_emax - 1.0, "KPM e_max too low");
}
#[test]
fn b2_kpm_dos_peak_locations() {
let (elems, pos) = ethanol();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let config = KpmConfig {
order: 200,
..Default::default()
};
let dos = compute_kpm_dos(&h, &config, -30.0, 10.0, 500);
let max_dos = dos
.total_dos
.iter()
.cloned()
.fold(f64::NEG_INFINITY, f64::max);
let peak_idx = dos.total_dos.iter().position(|&v| v == max_dos).unwrap();
let peak_energy = dos.energies[peak_idx];
let eht = solve_eht(&elems, &pos, None).unwrap();
let min_dist: f64 = eht
.energies
.iter()
.map(|&e| (e - peak_energy).abs())
.fold(f64::INFINITY, f64::min);
println!(
"[B2] KPM DOS peak at {:.2} eV, nearest eigenvalue {:.2} eV away",
peak_energy, min_dist
);
assert!(min_dist < 5.0, "KPM DOS peak too far from any eigenvalue");
}
#[test]
fn b2_kpm_mulliken_vs_eht_population() {
let (elems, pos) = water();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let eht = solve_eht(&elems, &pos, None).unwrap();
let nuclear_charges: Vec<f64> = elems.iter().map(|&z| z as f64).collect();
let config = KpmConfig {
order: 200,
..Default::default()
};
let kpm_mull = compute_kpm_mulliken(&h, &s, eht.n_electrons, &nuclear_charges, &config);
let pop = sci_form::compute_population(&elems, &pos).unwrap();
println!("[B2] EHT Mulliken: {:?}", pop.mulliken_charges);
println!("[B2] KPM Mulliken: {:?}", kpm_mull.charges);
let o_sign_agree =
kpm_mull.charges[0] * pop.mulliken_charges[0] >= 0.0 || kpm_mull.charges[0].abs() < 0.1;
println!("[B2] O charge sign agreement: {}", o_sign_agree);
}
#[test]
fn b2_timing_kpm_vs_exact_diag() {
let (elems, pos) = ethanol();
let basis = sci_form::eht::basis::build_basis(&elems, &pos);
let s = sci_form::eht::overlap::build_overlap_matrix(&basis);
let h = sci_form::eht::hamiltonian::build_hamiltonian(&basis, &s, None);
let t0 = Instant::now();
for _ in 0..50 {
let _ = solve_eht(&elems, &pos, None);
}
let eht_us = t0.elapsed().as_micros() as f64 / 50.0;
let config = KpmConfig {
order: 100,
..Default::default()
};
let t0 = Instant::now();
for _ in 0..50 {
let _ = compute_kpm_dos(&h, &config, -30.0, 10.0, 200);
}
let kpm_us = t0.elapsed().as_micros() as f64 / 50.0;
println!(
"[B2] Timing: EHT exact={:.0} µs, KPM DOS={:.0} µs",
eht_us, kpm_us
);
}
}
#[cfg(feature = "experimental-eeq")]
mod bench_eeq_vs_gasteiger {
use sci_form::charges_eeq::*;
use std::time::Instant;
fn water() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![8, 1, 1],
vec![
[0.0, 0.0, 0.117],
[0.0, 0.757, -0.469],
[0.0, -0.757, -0.469],
],
)
}
fn acetic_acid() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![6, 6, 8, 8, 1, 1, 1, 1],
vec![
[0.0, 0.0, 0.0],
[1.52, 0.0, 0.0],
[2.08, 1.20, 0.0],
[2.08, -1.20, 0.0],
[-0.36, 1.03, 0.0],
[-0.36, -0.51, 0.89],
[-0.36, -0.51, -0.89],
[3.05, -1.20, 0.0],
],
)
}
const WATER_REF_O: f64 = -0.80;
#[test]
fn b3_charge_sign_agreement() {
let (elems, pos) = water();
let gasteiger = sci_form::compute_charges("O").unwrap();
let eeq = compute_eeq_charges(&elems, &pos, &EeqConfig::default());
println!(
"[B3] Gasteiger O={:.4}, H={:.4}",
gasteiger.charges[0], gasteiger.charges[1]
);
println!(
"[B3] EEQ O={:.4}, H={:.4}",
eeq.charges[0], eeq.charges[1]
);
assert!(gasteiger.charges[0] < 0.0, "Gasteiger O should be negative");
assert!(eeq.charges[0] < 0.0, "EEQ O should be negative");
}
#[test]
fn b3_deviation_from_reference_water() {
let (elems, pos) = water();
let eeq = compute_eeq_charges(&elems, &pos, &EeqConfig::default());
let gasteiger = sci_form::compute_charges("O").unwrap();
let eeq_dev_o = (eeq.charges[0] - WATER_REF_O).abs();
let gas_dev_o = (gasteiger.charges[0] - WATER_REF_O).abs();
println!(
"[B3] Deviation from HF ref (O): EEQ={:.4}, Gasteiger={:.4}",
eeq_dev_o, gas_dev_o
);
assert!(eeq_dev_o < 1.5, "EEQ deviation too large");
assert!(gas_dev_o < 1.5, "Gasteiger deviation too large");
}
#[test]
fn b3_charge_neutrality() {
let (elems, pos) = acetic_acid();
let eeq = compute_eeq_charges(&elems, &pos, &EeqConfig::default());
let total: f64 = eeq.charges.iter().sum();
println!("[B3] EEQ total charge (acetic acid): {:.6}", total);
assert!(
total.abs() < 0.01,
"EEQ should preserve charge neutrality: {:.6}",
total
);
}
#[test]
fn b3_geometry_dependence() {
let (elems, _) = water();
let pos_normal = vec![
[0.0, 0.0, 0.117],
[0.0, 0.757, -0.469],
[0.0, -0.757, -0.469],
];
let pos_stretched = vec![
[0.0, 0.0, 0.117],
[0.0, 1.50, -0.469],
[0.0, -0.757, -0.469],
];
let eeq_normal = compute_eeq_charges(&elems, &pos_normal, &EeqConfig::default());
let eeq_stretched = compute_eeq_charges(&elems, &pos_stretched, &EeqConfig::default());
let delta_o = (eeq_normal.charges[0] - eeq_stretched.charges[0]).abs();
println!("[B3] EEQ O charge change on stretching: {:.4}", delta_o);
assert!(delta_o > 1e-4, "EEQ should be geometry-dependent");
}
#[test]
fn b3_timing_comparison() {
let (elems, pos) = acetic_acid();
let t0 = Instant::now();
for _ in 0..200 {
let _ = sci_form::compute_charges("CC(=O)O");
}
let gas_us = t0.elapsed().as_micros() as f64 / 200.0;
let config = EeqConfig::default();
let t0 = Instant::now();
for _ in 0..200 {
let _ = compute_eeq_charges(&elems, &pos, &config);
}
let eeq_us = t0.elapsed().as_micros() as f64 / 200.0;
println!(
"[B3] Timing: Gasteiger={:.0} µs, EEQ={:.0} µs",
gas_us, eeq_us
);
}
}
#[cfg(feature = "experimental-alpb")]
mod bench_alpb_vs_gb {
use sci_form::solvation::{compute_gb_solvation, compute_nonpolar_solvation};
use sci_form::solvation_alpb::*;
use std::time::Instant;
fn water() -> (Vec<u8>, Vec<[f64; 3]>, Vec<f64>) {
(
vec![8, 1, 1],
vec![
[0.0, 0.0, 0.117],
[0.0, 0.757, -0.469],
[0.0, -0.757, -0.469],
],
vec![-0.834, 0.417, 0.417],
)
}
fn ethanol() -> (Vec<u8>, Vec<[f64; 3]>, Vec<f64>) {
(
vec![6, 6, 8, 1, 1, 1, 1, 1, 1],
vec![
[0.0, 0.0, 0.0],
[1.52, 0.0, 0.0],
[2.13, 1.21, 0.0],
[-0.36, 1.03, 0.0],
[-0.36, -0.51, 0.89],
[-0.36, -0.51, -0.89],
[1.88, -0.51, 0.89],
[1.88, -0.51, -0.89],
[3.09, 1.21, 0.0],
],
vec![-0.10, 0.15, -0.68, 0.04, 0.04, 0.04, 0.06, 0.06, 0.39],
)
}
const WATER_SOLVATION_REF: f64 = -6.3;
#[test]
fn b4_solvation_sign_agreement() {
let (elems, pos, charges) = water();
let gb = compute_gb_solvation(&elems, &pos, &charges, Some(78.5), Some(1.0), Some(1.4));
let alpb = compute_alpb_solvation(&elems, &pos, &charges, &AlpbConfig::default());
println!(
"[B4] GB-HCT total: {:.3} kcal/mol",
gb.total_energy_kcal_mol
);
println!("[B4] ALPB total: {:.3} kcal/mol", alpb.total_energy);
assert!(
gb.total_energy_kcal_mol < 0.0 || gb.electrostatic_energy_kcal_mol < 0.0,
"GB should give negative solvation"
);
assert!(
alpb.total_energy < 0.0 || alpb.electrostatic_energy < 0.0,
"ALPB should give negative solvation"
);
}
#[test]
fn b4_deviation_from_experimental() {
let (elems, pos, charges) = water();
let gb = compute_gb_solvation(&elems, &pos, &charges, Some(78.5), Some(1.0), Some(1.4));
let alpb = compute_alpb_solvation(&elems, &pos, &charges, &AlpbConfig::default());
let gb_dev = (gb.total_energy_kcal_mol - WATER_SOLVATION_REF).abs();
let alpb_dev = (alpb.total_energy - WATER_SOLVATION_REF).abs();
println!(
"[B4] Deviation from exp ({:.1} kcal/mol): GB={:.2}, ALPB={:.2}",
WATER_SOLVATION_REF, gb_dev, alpb_dev
);
assert!(gb_dev < 30.0, "GB deviation too large");
assert!(alpb_dev < 30.0, "ALPB deviation too large");
}
#[test]
fn b4_born_radii_comparison() {
let (elems, pos, _) = water();
let gb_born = sci_form::solvation::compute_born_radii(&elems, &pos);
let alpb_born = compute_born_radii(&elems, &pos, 1.4);
println!("[B4] GB-HCT Born radii: {:?}", gb_born);
println!("[B4] ALPB Born radii: {:?}", alpb_born.radii);
for (i, (&gb_r, &alpb_r)) in gb_born.iter().zip(alpb_born.radii.iter()).enumerate() {
assert!(gb_r > 0.0, "GB Born[{}]={}", i, gb_r);
assert!(alpb_r > 0.3 && alpb_r < 50.0, "ALPB Born[{}]={}", i, alpb_r);
}
}
#[test]
fn b4_nonpolar_component_agreement() {
let (elems, pos, _) = water();
let np = compute_nonpolar_solvation(&elems, &pos, None);
let charges = vec![-0.834, 0.417, 0.417];
let alpb = compute_alpb_solvation(&elems, &pos, &charges, &AlpbConfig::default());
println!("[B4] Legacy non-polar: {:.4} kcal/mol", np.energy_kcal_mol);
println!(
"[B4] ALPB non-polar: {:.4} kcal/mol",
alpb.nonpolar_energy
);
assert!(np.energy_kcal_mol.abs() < 20.0, "Legacy NP too large");
assert!(alpb.nonpolar_energy.abs() < 20.0, "ALPB NP too large");
}
#[test]
fn b4_timing_comparison() {
let (elems, pos, charges) = ethanol();
let t0 = Instant::now();
for _ in 0..100 {
let _ = compute_gb_solvation(&elems, &pos, &charges, Some(78.5), Some(1.0), Some(1.4));
}
let gb_us = t0.elapsed().as_micros() as f64 / 100.0;
let config = AlpbConfig::default();
let t0 = Instant::now();
for _ in 0..100 {
let _ = compute_alpb_solvation(&elems, &pos, &charges, &config);
}
let alpb_us = t0.elapsed().as_micros() as f64 / 100.0;
println!(
"[B4] Timing: GB-HCT={:.0} µs, ALPB={:.0} µs",
gb_us, alpb_us
);
}
}
#[cfg(feature = "experimental-d4")]
mod bench_d4_vs_uff {
use sci_form::dispersion::*;
use std::time::Instant;
fn benzene_dimer() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![
6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1, 6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1,
],
vec![
[1.40, 0.0, 0.0],
[0.70, 1.21, 0.0],
[-0.70, 1.21, 0.0],
[-1.40, 0.0, 0.0],
[-0.70, -1.21, 0.0],
[0.70, -1.21, 0.0],
[2.48, 0.0, 0.0],
[1.24, 2.15, 0.0],
[-1.24, 2.15, 0.0],
[-2.48, 0.0, 0.0],
[-1.24, -2.15, 0.0],
[1.24, -2.15, 0.0],
[1.40, 0.0, 3.5],
[0.70, 1.21, 3.5],
[-0.70, 1.21, 3.5],
[-1.40, 0.0, 3.5],
[-0.70, -1.21, 3.5],
[0.70, -1.21, 3.5],
[2.48, 0.0, 3.5],
[1.24, 2.15, 3.5],
[-1.24, 2.15, 3.5],
[-2.48, 0.0, 3.5],
[-1.24, -2.15, 3.5],
[1.24, -2.15, 3.5],
],
)
}
fn ethanol() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![6, 6, 8, 1, 1, 1, 1, 1, 1],
vec![
[0.0, 0.0, 0.0],
[1.52, 0.0, 0.0],
[2.13, 1.21, 0.0],
[-0.36, 1.03, 0.0],
[-0.36, -0.51, 0.89],
[-0.36, -0.51, -0.89],
[1.88, -0.51, 0.89],
[1.88, -0.51, -0.89],
[3.09, 1.21, 0.0],
],
)
}
const BENZENE_DIMER_REF_KCAL: f64 = -2.7;
#[test]
fn b5_d4_dispersion_attractive() {
let (elems, pos) = benzene_dimer();
let d4 = compute_d4_energy(&elems, &pos, &D4Config::default());
println!(
"[B5] D4 total: {:.4} Hartree = {:.2} kcal/mol",
d4.total_energy, d4.total_kcal_mol
);
assert!(
d4.total_energy < 0.0,
"D4 dispersion should be negative (attractive)"
);
}
#[test]
fn b5_d4_vs_reference_benzene_dimer() {
let (elems, pos) = benzene_dimer();
let d4 = compute_d4_energy(&elems, &pos, &D4Config::default());
println!(
"[B5] D4 dimer: {:.2} kcal/mol, ref: {:.2} kcal/mol, Δ={:.2}",
d4.total_kcal_mol,
BENZENE_DIMER_REF_KCAL,
(d4.total_kcal_mol - BENZENE_DIMER_REF_KCAL).abs()
);
assert!(d4.total_kcal_mol < 0.0, "D4 should be attractive");
assert!(d4.total_kcal_mol.abs() < 50.0, "D4 magnitude reasonable");
}
#[test]
fn b5_uff_augmented_with_d4() {
let (elems, pos) = ethanol();
let conf = sci_form::embed("CCO", 42);
let uff_energy = sci_form::compute_uff_energy("CCO", &conf.coords).unwrap_or(0.0);
let d4 = compute_d4_energy(&elems, &pos, &D4Config::default());
let augmented = uff_energy + d4.total_kcal_mol;
println!(
"[B5] UFF: {:.2}, D4: {:.4}, UFF+D4: {:.2} kcal/mol",
uff_energy, d4.total_kcal_mol, augmented
);
assert!(augmented <= uff_energy + 0.01, "D4 should lower UFF energy");
}
#[test]
fn b5_coordination_numbers() {
let (elems, pos) = ethanol();
let d4 = compute_d4_energy(&elems, &pos, &D4Config::default());
println!(
"[B5] D4 coordination numbers: {:?}",
d4.coordination_numbers
);
assert!(d4.coordination_numbers[0] > 2.0, "C CN too low");
assert!(d4.coordination_numbers[2] > 0.5, "O CN too low");
}
#[test]
fn b5_timing_d4_overhead() {
let (elems, pos) = ethanol();
let conf = sci_form::embed("CCO", 42);
let t0 = Instant::now();
for _ in 0..200 {
let _ = sci_form::compute_uff_energy("CCO", &conf.coords);
}
let uff_us = t0.elapsed().as_micros() as f64 / 200.0;
let t0 = Instant::now();
for _ in 0..200 {
let _ = compute_d4_energy(&elems, &pos, &D4Config::default());
}
let d4_us = t0.elapsed().as_micros() as f64 / 200.0;
println!(
"[B5] Timing: UFF={:.0} µs, D4={:.0} µs, overhead={:.1}%",
uff_us,
d4_us,
d4_us / uff_us.max(1.0) * 100.0
);
}
}
#[cfg(feature = "experimental-riemannian")]
mod bench_riemannian_vs_etkdg {
use nalgebra::DMatrix;
use sci_form::beta::riemannian::*;
use std::time::Instant;
fn make_constraints(coords: &[[f64; 3]]) -> Vec<DistanceConstraint> {
let n = coords.len();
let mut constraints = Vec::new();
for i in 0..n {
for j in (i + 1)..n {
let dx = coords[i][0] - coords[j][0];
let dy = coords[i][1] - coords[j][1];
let dz = coords[i][2] - coords[j][2];
let d = (dx * dx + dy * dy + dz * dz).sqrt();
constraints.push(DistanceConstraint {
i,
j,
lower: d * 0.95,
upper: d * 1.05,
});
}
}
constraints
}
#[test]
fn b6_psd_guarantee() {
let coords = vec![
[0.0, 0.0, 0.0],
[1.5, 0.0, 0.0],
[0.75, 1.3, 0.0],
[0.75, 0.43, 1.22],
];
let n = coords.len();
let constraints = make_constraints(&coords);
let config = RiemannianConfig::default();
let lbfgs = RiemannianLbfgs::new(config);
let initial = DMatrix::identity(n, n);
let result = lbfgs.optimize(&initial, &constraints);
println!(
"[B6] Riemannian converged: {}, iter: {}, max_violation: {:.6}",
result.converged, result.iterations, result.max_violation
);
let eig = nalgebra::SymmetricEigen::new(result.metric_matrix.clone());
let min_eig = eig
.eigenvalues
.iter()
.cloned()
.fold(f64::INFINITY, f64::min);
println!("[B6] Min eigenvalue of metric matrix: {:.6e}", min_eig);
assert!(
min_eig >= -1e-10,
"Metric matrix not PSD: min eigenvalue = {}",
min_eig
);
}
#[test]
fn b6_etkdg_vs_riemannian_simple() {
let conf = sci_form::embed("C", 42);
assert!(conf.error.is_none(), "embed failed: {:?}", conf.error);
let n = conf.num_atoms;
let pos: Vec<[f64; 3]> = conf.coords.chunks(3).map(|c| [c[0], c[1], c[2]]).collect();
let constraints = make_constraints(&pos);
let config = RiemannianConfig::default();
let lbfgs = RiemannianLbfgs::new(config);
let initial = DMatrix::identity(n, n);
let result = lbfgs.optimize(&initial, &constraints);
println!("[B6] ETKDG time: {:.2} ms", conf.time_ms);
println!(
"[B6] Riemannian: converged={}, violation={:.4} Å, iter={}",
result.converged, result.max_violation, result.iterations
);
}
#[test]
fn b6_timing_comparison() {
let t0 = Instant::now();
for _ in 0..20 {
let _ = sci_form::embed("CCO", 42);
}
let etkdg_us = t0.elapsed().as_micros() as f64 / 20.0;
let n = 4;
let constraints = vec![
DistanceConstraint {
i: 0,
j: 1,
lower: 1.4,
upper: 1.6,
},
DistanceConstraint {
i: 0,
j: 2,
lower: 1.4,
upper: 1.6,
},
DistanceConstraint {
i: 0,
j: 3,
lower: 1.4,
upper: 1.6,
},
DistanceConstraint {
i: 1,
j: 2,
lower: 1.4,
upper: 1.6,
},
DistanceConstraint {
i: 1,
j: 3,
lower: 1.4,
upper: 1.6,
},
DistanceConstraint {
i: 2,
j: 3,
lower: 1.4,
upper: 1.6,
},
];
let config = RiemannianConfig::default();
let t0 = Instant::now();
for _ in 0..20 {
let lbfgs = RiemannianLbfgs::new(config.clone());
let initial = DMatrix::identity(n, n);
let _ = lbfgs.optimize(&initial, &constraints);
}
let riem_us = t0.elapsed().as_micros() as f64 / 20.0;
println!(
"[B6] Timing: ETKDG={:.0} µs, Riemannian={:.0} µs",
etkdg_us, riem_us
);
}
}
#[cfg(feature = "experimental-sdr")]
mod bench_sdr_vs_etkdg {
use sci_form::alpha::sdr::*;
use std::time::Instant;
#[test]
fn b7_sdr_embedding_basic() {
let n = 5;
let d_pairs = vec![
(0, 1, 1.5),
(1, 2, 1.5),
(2, 3, 1.5),
(3, 4, 1.5),
(0, 2, 2.45),
(1, 3, 2.45),
(2, 4, 2.45),
];
let result = sdr_embed(n, &d_pairs, &SdrConfig::default());
println!(
"[B7] SDR max_distance_error: {:.4} Å, retries_avoided: {}",
result.max_distance_error, result.retries_avoided
);
assert_eq!(result.num_atoms, n);
assert_eq!(result.coords.len(), n * 3);
}
#[test]
fn b7_sdr_vs_etkdg_distance_fidelity() {
let conf = sci_form::embed("CC", 42);
assert!(conf.error.is_none());
let n = conf.num_atoms;
let pos: Vec<[f64; 3]> = conf.coords.chunks(3).map(|c| [c[0], c[1], c[2]]).collect();
let mut d_pairs = Vec::new();
for i in 0..n {
for j in (i + 1)..n {
let dx = pos[i][0] - pos[j][0];
let dy = pos[i][1] - pos[j][1];
let dz = pos[i][2] - pos[j][2];
let d = (dx * dx + dy * dy + dz * dz).sqrt();
d_pairs.push((i, j, d));
}
}
let sdr = sdr_embed(n, &d_pairs, &SdrConfig::default());
let sdr_pos: Vec<[f64; 3]> = sdr.coords.chunks(3).map(|c| [c[0], c[1], c[2]]).collect();
let mut max_err = 0.0_f64;
for &(i, j, d_target) in &d_pairs {
let dx = sdr_pos[i][0] - sdr_pos[j][0];
let dy = sdr_pos[i][1] - sdr_pos[j][1];
let dz = sdr_pos[i][2] - sdr_pos[j][2];
let d_sdr = (dx * dx + dy * dy + dz * dz).sqrt();
max_err = max_err.max((d_sdr - d_target).abs());
}
println!("[B7] ETKDG → SDR max distance error: {:.4} Å", max_err);
assert!(
max_err < 1.0,
"SDR distance fidelity too poor: {:.3} Å",
max_err
);
}
#[test]
fn b7_timing_comparison() {
let t0 = Instant::now();
for _ in 0..20 {
let _ = sci_form::embed("CCO", 42);
}
let etkdg_us = t0.elapsed().as_micros() as f64 / 20.0;
let n = 5;
let d_pairs = vec![
(0, 1, 1.5),
(1, 2, 1.5),
(2, 3, 1.5),
(3, 4, 1.5),
(0, 2, 2.45),
(1, 3, 2.45),
(2, 4, 2.45),
];
let t0 = Instant::now();
for _ in 0..20 {
let _ = sdr_embed(n, &d_pairs, &SdrConfig::default());
}
let sdr_us = t0.elapsed().as_micros() as f64 / 20.0;
println!(
"[B7] Timing: ETKDG={:.0} µs, SDR={:.0} µs",
etkdg_us, sdr_us
);
}
}
#[cfg(feature = "experimental-mbh")]
mod bench_mbh_vs_full_hessian {
use sci_form::beta::mbh::*;
use std::time::Instant;
fn ethanol() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![6, 6, 8, 1, 1, 1, 1, 1, 1],
vec![
[0.0, 0.0, 0.0],
[1.52, 0.0, 0.0],
[2.13, 1.21, 0.0],
[-0.36, 1.03, 0.0],
[-0.36, -0.51, 0.89],
[-0.36, -0.51, -0.89],
[1.88, -0.51, 0.89],
[1.88, -0.51, -0.89],
[3.09, 1.21, 0.0],
],
)
}
fn benzene() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1],
vec![
[1.40, 0.0, 0.0],
[0.70, 1.21, 0.0],
[-0.70, 1.21, 0.0],
[-1.40, 0.0, 0.0],
[-0.70, -1.21, 0.0],
[0.70, -1.21, 0.0],
[2.48, 0.0, 0.0],
[1.24, 2.15, 0.0],
[-1.24, 2.15, 0.0],
[-2.48, 0.0, 0.0],
[-1.24, -2.15, 0.0],
[1.24, -2.15, 0.0],
],
)
}
fn uff_energy_fn(smiles: &'static str) -> impl Fn(&[f64]) -> f64 {
move |coords: &[f64]| -> f64 { sci_form::compute_uff_energy(smiles, coords).unwrap_or(0.0) }
}
#[test]
fn b8_mbh_dof_reduction() {
let (elems, pos) = benzene();
let rings: Vec<(Vec<usize>, bool)> = vec![(vec![0, 1, 2, 3, 4, 5], true)];
let energy_fn = uff_energy_fn("c1ccccc1");
let mbh = compute_mbh_frequencies(&elems, &pos, &rings, &energy_fn, &MbhConfig::default());
println!(
"[B8] Full DOF: {}, Reduced DOF: {}, Speedup: {:.1}x",
mbh.n_dof_full, mbh.n_dof_reduced, mbh.speedup
);
assert!(mbh.speedup >= 1.0, "Speedup should be >= 1.0");
assert!(mbh.n_dof_reduced <= mbh.n_dof_full);
}
#[test]
fn b8_mbh_frequency_physical() {
let (elems, pos) = ethanol();
let rings: Vec<(Vec<usize>, bool)> = vec![];
let energy_fn = uff_energy_fn("CCO");
let mbh = compute_mbh_frequencies(&elems, &pos, &rings, &energy_fn, &MbhConfig::default());
println!(
"[B8] MBH frequencies (N={}): {:?}",
mbh.frequencies.len(),
mbh.frequencies
);
let positive: Vec<_> = mbh.frequencies.iter().filter(|&&f| f > 10.0).collect();
assert!(
!positive.is_empty(),
"Should have positive vibrational frequencies"
);
}
#[test]
fn b8_mbh_vs_full_ethanol() {
let (elems, pos) = ethanol();
let full = sci_form::ir::vibrations::compute_vibrational_analysis(
&elems,
&pos,
sci_form::ir::hessian::HessianMethod::Pm3,
None,
);
let rings: Vec<(Vec<usize>, bool)> = vec![];
let energy_fn = uff_energy_fn("CCO");
let mbh = compute_mbh_frequencies(&elems, &pos, &rings, &energy_fn, &MbhConfig::default());
if let Ok(ref full_result) = full {
let full_freqs: Vec<f64> = full_result
.modes
.iter()
.filter(|m| m.is_real && m.frequency_cm1 > 10.0)
.map(|m| m.frequency_cm1)
.collect();
println!(
"[B8] Full Hessian (N={}): [{:.0} ... {:.0}]",
full_freqs.len(),
full_freqs.first().unwrap_or(&0.0),
full_freqs.last().unwrap_or(&0.0)
);
}
println!(
"[B8] MBH (N={}): [{:.0} ... {:.0}]",
mbh.frequencies.len(),
mbh.frequencies.first().unwrap_or(&0.0),
mbh.frequencies.last().unwrap_or(&0.0)
);
}
#[test]
fn b8_timing_comparison() {
let (elems, pos) = ethanol();
let rings: Vec<(Vec<usize>, bool)> = vec![];
let energy_fn = uff_energy_fn("CCO");
let t0 = Instant::now();
for _ in 0..3 {
let _ =
compute_mbh_frequencies(&elems, &pos, &rings, &energy_fn, &MbhConfig::default());
}
let mbh_us = t0.elapsed().as_micros() as f64 / 3.0;
let t0 = Instant::now();
let _ = sci_form::ir::vibrations::compute_vibrational_analysis(
&elems,
&pos,
sci_form::ir::hessian::HessianMethod::Pm3,
None,
);
let full_us = t0.elapsed().as_micros() as f64;
println!(
"[B8] Timing: Full Hessian={:.0} µs, MBH={:.0} µs",
full_us, mbh_us
);
}
}
#[cfg(feature = "experimental-cga")]
mod bench_cga_geometry {
use sci_form::alpha::cga::*;
#[test]
fn b9_motor_rotation_preserves_distance() {
let a = [0.0, 0.0, 0.0];
let b = [1.5, 0.0, 0.0];
let angle = std::f64::consts::PI / 4.0;
let motor = dihedral_motor(a, b, angle);
let p = [2.0, 1.0, 0.0];
let p_rot = motor.transform_point(p);
let dist_orig = (p[1] * p[1] + p[2] * p[2]).sqrt();
let dist_rot = (p_rot[1] * p_rot[1] + p_rot[2] * p_rot[2]).sqrt();
println!(
"[B9] Original dist from axis: {:.4}, Rotated: {:.4}",
dist_orig, dist_rot
);
assert!(
(dist_orig - dist_rot).abs() < 0.1,
"Motor should preserve distance"
);
}
#[test]
fn b9_motor_identity() {
let motor = Motor::identity();
let p = [1.0, 2.0, 3.0];
let result = motor.transform_point(p);
assert!((result[0] - p[0]).abs() < 1e-10);
assert!((result[1] - p[1]).abs() < 1e-10);
assert!((result[2] - p[2]).abs() < 1e-10);
}
#[test]
fn b9_cga_vs_matrix_rotation() {
let axis_a = [0.0, 0.0, 0.0];
let axis_b = [0.0, 0.0, 1.0];
let angle = std::f64::consts::PI / 2.0;
let motor = dihedral_motor(axis_a, axis_b, angle);
let p = [1.0, 0.0, 0.0];
let result = motor.transform_point(p);
println!(
"[B9] Rotated [1,0,0] by 90° about z: [{:.3}, {:.3}, {:.3}]",
result[0], result[1], result[2]
);
assert!(result[0].abs() < 0.2, "x should be ~0");
assert!((result[1] - 1.0).abs() < 0.2, "y should be ~1");
assert!(result[2].abs() < 0.2, "z should be ~0");
}
#[test]
fn b9_subtree_rotation() {
let coords = vec![0.0, 0.0, 0.0, 1.5, 0.0, 0.0, 3.0, 0.0, 0.0];
let a = [0.0, 0.0, 0.0];
let b = [1.5, 0.0, 0.0];
let motor = dihedral_motor(a, b, std::f64::consts::PI / 2.0);
let result = apply_motor_to_subtree(&coords, &[2], &motor);
println!(
"[B9] After subtree rotation: atom2 = [{:.3}, {:.3}, {:.3}]",
result[6], result[7], result[8]
);
assert!((result[0] - coords[0]).abs() < 1e-10);
assert!((result[3] - coords[3]).abs() < 1e-10);
}
}
#[cfg(feature = "experimental-cpm")]
mod bench_cpm_electrochemistry {
use sci_form::beta::cpm::*;
fn water() -> (Vec<u8>, Vec<[f64; 3]>) {
(
vec![8, 1, 1],
vec![
[0.0, 0.0, 0.117],
[0.0, 0.757, -0.469],
[0.0, -0.757, -0.469],
],
)
}
#[test]
fn b10_cpm_charge_varies_with_potential() {
let (elems, pos) = water();
let mut charges_at_potentials = Vec::new();
for mu in [-5.0_f64, -4.5, -4.0, -3.5] {
let config = CpmConfig {
mu_ev: mu,
..Default::default()
};
let result = compute_cpm_charges(&elems, &pos, &config);
charges_at_potentials.push((mu, result.total_charge));
println!("[B10] µ = {:.1} eV → Q = {:.4} e", mu, result.total_charge);
}
let q_max = charges_at_potentials
.iter()
.map(|&(_, q)| q)
.fold(f64::NEG_INFINITY, f64::max);
let q_min = charges_at_potentials
.iter()
.map(|&(_, q)| q)
.fold(f64::INFINITY, f64::min);
let q_range = q_max - q_min;
println!("[B10] Charge range: {:.4} e", q_range);
assert!(q_range > 1e-6, "Charge should vary with potential");
}
#[test]
fn b10_cpm_surface_shape() {
let (elems, pos) = water();
let surface = compute_cpm_surface(&elems, &pos, -5.5, -3.5, 10, 78.5);
println!("[B10] CPM surface points: {}", surface.mu_values.len());
assert_eq!(surface.mu_values.len(), 10);
assert_eq!(surface.total_charge.len(), 10);
assert_eq!(surface.free_energy.len(), 10);
}
#[test]
fn b10_grand_potential_decreases() {
let (elems, pos) = water();
let surface = compute_cpm_surface(&elems, &pos, -5.5, -3.5, 20, 78.5);
let finite_energies: Vec<f64> = surface
.free_energy
.iter()
.filter(|e| e.is_finite())
.copied()
.collect();
if finite_energies.len() >= 2 {
println!(
"[B10] Free energy range: [{:.4}, {:.4}] eV",
finite_energies
.iter()
.cloned()
.fold(f64::INFINITY, f64::min),
finite_energies
.iter()
.cloned()
.fold(f64::NEG_INFINITY, f64::max)
);
}
assert!(
surface.all_converged,
"CPM surface should converge at all points"
);
}
#[test]
fn b10_capacitance_positive() {
let (elems, pos) = water();
let surface = compute_cpm_surface(&elems, &pos, -5.5, -3.5, 10, 78.5);
for (i, &c) in surface.capacitance.iter().enumerate() {
println!("[B10] C(µ={:.2}) = {:.4} e/eV", surface.mu_values[i], c);
}
}
}
#[cfg(feature = "experimental-gsm")]
mod bench_gsm_reaction_path {
use sci_form::alpha::gsm::*;
fn h2_energy(coords: &[f64]) -> f64 {
let dx = coords[3] - coords[0];
let dy = coords[4] - coords[1];
let dz = coords[5] - coords[2];
let d = (dx * dx + dy * dy + dz * dz).sqrt();
let d0 = 0.74;
let k = 200.0;
k * (d - d0).powi(2)
}
#[test]
fn b11_string_growth() {
let reactant = vec![0.0, 0.0, 0.0, 0.74, 0.0, 0.0];
let product = vec![0.0, 0.0, 0.0, 3.0, 0.0, 0.0];
let config = GsmConfig {
max_nodes: 7,
max_iter: 20,
..Default::default()
};
let path = grow_string(&reactant, &product, &h2_energy, &config);
println!(
"[B11] GSM nodes: {}, joined: {}",
path.nodes.len(),
path.joined
);
println!("[B11] Energies: {:?}", path.energies);
assert!(
path.nodes.len() >= 2,
"Should have at least reactant and product"
);
}
#[test]
fn b11_interpolation_smoothness() {
let reactant = vec![0.0, 0.0, 0.0, 0.74, 0.0, 0.0];
let product = vec![0.0, 0.0, 0.0, 3.0, 0.0, 0.0];
let mut prev_d = 0.0_f64;
for &t in &[0.0, 0.25, 0.5, 0.75, 1.0] {
let node = interpolate_node(&reactant, &product, t);
let d = ((node[3] - node[0]).powi(2)
+ (node[4] - node[1]).powi(2)
+ (node[5] - node[2]).powi(2))
.sqrt();
println!("[B11] t={:.2} → d(H-H) = {:.3} Å", t, d);
if t > 0.0 {
assert!(d >= prev_d - 0.01, "Distance should increase monotonically");
}
prev_d = d;
}
}
#[test]
fn b11_transition_state_search() {
let reactant = vec![0.0, 0.0, 0.0, 0.74, 0.0, 0.0];
let product = vec![0.0, 0.0, 0.0, 3.0, 0.0, 0.0];
let config = GsmConfig {
max_nodes: 7,
max_iter: 10,
..Default::default()
};
let result = find_transition_state(&reactant, &product, &h2_energy, &config);
println!(
"[B11] TS energy: {:.2} kcal/mol, activation: {:.2} kcal/mol",
result.ts_energy, result.activation_energy
);
println!(
"[B11] TS node index: {}, n_nodes: {}",
result.ts_node_index, result.n_nodes
);
assert!(
result.activation_energy >= 0.0,
"Activation energy should be non-negative"
);
}
#[test]
fn b11_energy_profile_has_maximum() {
let reactant = vec![0.0, 0.0, 0.0, 0.74, 0.0, 0.0];
let product = vec![0.0, 0.0, 0.0, 3.0, 0.0, 0.0];
let config = GsmConfig {
max_nodes: 11,
max_iter: 20,
..Default::default()
};
let result = find_transition_state(&reactant, &product, &h2_energy, &config);
if result.path_energies.len() >= 3 {
let e_first = result.path_energies[0];
let e_last = *result.path_energies.last().unwrap();
let e_max = result
.path_energies
.iter()
.cloned()
.fold(f64::NEG_INFINITY, f64::max);
println!(
"[B11] E(react)={:.2}, E(max)={:.2}, E(prod)={:.2}",
e_first, e_max, e_last
);
assert!(
e_max >= e_first - 0.01,
"Maximum should be >= reactant energy"
);
}
}
}