interatomic 0.6.1

//! Investigate how few spline coefficients are needed for SR potentials.
//!
//! Run with: `cargo run --example spline_reduction --release`

use interatomic::twobody::*;
use interatomic::Cutoff;

/// Max absolute errors over 100k test points in [r_min_test, cutoff].
fn validate_abs<P: IsotropicTwobodyEnergy + Cutoff>(
    potential: &P,
    splined: &SplinedPotential,
    r_min_test: f64,
) -> (f64, f64) {
    let rsq_min = r_min_test * r_min_test;
    let rsq_max = splined.cutoff_squared();
    let n = 100_000;
    let (mut max_u, mut max_f) = (0.0f64, 0.0f64);
    let shift = splined.stats().energy_shift;
    for i in 0..n {
        let t = (i as f64 + 0.37) / n as f64;
        let rsq = rsq_min + t * (rsq_max - rsq_min);
        max_u = max_u.max(
            (splined.isotropic_twobody_energy(rsq)
                - (potential.isotropic_twobody_energy(rsq) - shift))
                .abs(),
        );
        max_f = max_f.max(
            (splined.isotropic_twobody_force(rsq) - potential.isotropic_twobody_force(rsq)).abs(),
        );
    }
    (max_u, max_f)
}

fn analyze<P: IsotropicTwobodyEnergy + Cutoff>(
    name: &str,
    potential: &P,
    cutoff: f64,
    r_min_test: f64,
) {
    let grid_sizes = [50, 75, 100, 150, 200, 300, 500];
    println!("\n{name} (r > {r_min_test:.1} Å)");
    println!(
        "{:<8} {:>14} {:>14} {:>8}",
        "n_pts", "|dU| kJ/mol", "|dF| kJ/Å", "KB(f32)"
    );
    for &n in &grid_sizes {
        let config = SplineConfig {
            n_points: n,
            grid_type: GridType::PowerLaw2,
            shift_energy: true,
            ..Default::default()
        };
        let splined = SplinedPotential::with_cutoff(potential, cutoff, config);
        let (max_u, max_f) = validate_abs(potential, &splined, r_min_test);
        // GPU/SIMD: 8 x f32 per interval = 32 bytes
        let kb_f32 = n as f64 * 32.0 / 1024.0;
        println!("{:<8} {:>14.4e} {:>14.4e} {:>8.1}", n, max_u, max_f, kb_f32);
    }
}

fn main() {
    let cutoff = 20.0;

    // Smallest sigma in CALVADOS/KH is ~4 Å. Test from there.
    let sigma_min = 4.0;

    // === Ashbaugh-Hatch (CALVADOS) ===
    let lj1 = LennardJones::new(0.8368, 5.04);
    let ah = AshbaughHatch::new(lj1, 0.5, cutoff);
    analyze(
        "AshbaughHatch (ε=0.84, σ=5.04, λ=0.5)",
        &ah,
        cutoff,
        sigma_min,
    );

    let lj2 = LennardJones::new(0.8368, 5.04);
    let ah_rep = AshbaughHatch::new(lj2, 0.0, cutoff);
    analyze(
        "AshbaughHatch WCA (ε=0.84, σ=5.04, λ=0)",
        &ah_rep,
        cutoff,
        sigma_min,
    );

    // Smallest sigma case
    let lj_small = LennardJones::new(0.8368, 4.0);
    let ah_small = AshbaughHatch::new(lj_small, 0.5, cutoff);
    analyze(
        "AshbaughHatch (ε=0.84, σ=4.0, λ=0.5)",
        &ah_small,
        cutoff,
        sigma_min,
    );

    // === Kim-Hummer ===
    let kh_att = KimHummer::new(-0.5, 6.0);
    analyze(
        "KimHummer attractive (ε=-0.5, σ=6.0)",
        &kh_att,
        cutoff,
        sigma_min,
    );

    let kh_rep = KimHummer::new(0.5, 6.0);
    analyze(
        "KimHummer repulsive (ε=+0.5, σ=6.0)",
        &kh_rep,
        cutoff,
        sigma_min,
    );

    // KH with smallest sigma
    let kh_small = KimHummer::new(-0.5, 4.0);
    analyze("KimHummer (ε=-0.5, σ=4.0)", &kh_small, cutoff, sigma_min);

    // === LJ ===
    let lj = LennardJones::new(1.0, 4.0);
    analyze("LennardJones (ε=1.0, σ=4.0)", &lj, cutoff, sigma_min);

    // === Memory for typical systems ===
    println!("\n--- Memory per pair type (f32, GPU/SIMD) ---");
    println!("{:<8} {:>12} {:>12}", "n_pts", "per_pair", "400 pairs");
    for &n in &[100, 150, 200, 300, 500] {
        let per_pair = n * 32;
        let total = per_pair * 400;
        println!(
            "{:<8} {:>10.1} KB {:>10.1} KB",
            n,
            per_pair as f64 / 1024.0,
            total as f64 / 1024.0
        );
    }
}