chemx-props 0.3.0

//! Central-difference numerical Hessian from any gradient function.

use chemx_core::Molecule;
use rayon::prelude::*;

/// Compute the 3N×3N nuclear Hessian ∂²E/∂R² by central finite differences.
///
/// `gradient_fn(mol)` returns the flat gradient `[dE/dx₁, dE/dy₁, dE/dz₁, …,
/// dE/dz_N]` (length 3N). The function is called 2×3N times (one pair per
/// Cartesian degree of freedom). The returned Hessian is symmetrized
/// (H = (H + Hᵀ)/2) to suppress finite-difference noise.
///
/// **Step size.** 0.005 bohr gives ~10⁻¹⁰ Eh accuracy for RHF with tight SCF
/// (1×10⁻¹⁰ convergence); 0.01 bohr is faster with negligible extra error.
///
/// **Geometry.** Run at a geometry where the analytic gradient is ≤ ~10⁻⁵
/// hartree/bohr; modes will be spuriously imaginary at a stationary point that
/// has been under-converged.
///
/// **DFT grid noise.** Numerical Hessians of Kohn–Sham energies inherit the
/// quadrature noise of the XC grid: the gradient at each displaced geometry is
/// evaluated on that geometry's own (atom-following) grid, so low-frequency
/// modes can carry a few cm⁻¹ of grid noise at the default grid level. The
/// method's default grid is used deliberately (no silent grid bump); raise the
/// grid level if low modes look noisy.
pub fn numerical_hessian<F>(molecule: &Molecule, step_size: f64, gradient_fn: F) -> Vec<f64>
where
    F: Fn(&Molecule) -> Vec<f64> + Sync,
{
    let natom = molecule.len();
    let ndof = 3 * natom;

    // The 2·ndof displaced-geometry gradients are mutually independent, so they run
    // concurrently on the rayon pool (`gradient_fn` is `Fn + Sync` — each call builds
    // its own scratch / converges its own SCF). Each half-row of the raw Hessian comes
    // from exactly one (forward, backward) pair, assembled in a fixed order below, so
    // the result is bit-identical to the serial sweep — only the wall-clock changes.
    // This is the dominant cost of a frequency job (an N-atom Hessian is 6N gradients);
    // for small molecules the per-displacement gradients themselves barely use the
    // inner parallelism, so fanning the outer loop out is where the cores go to work.
    let evals: Vec<(usize, f64)> = (0..ndof)
        .flat_map(|dof| [(dof, step_size), (dof, -step_size)])
        .collect();
    let grads: Vec<Vec<f64>> = evals
        .par_iter()
        .map(|&(dof, delta)| gradient_fn(&displaced(molecule, dof / 3, dof % 3, delta)))
        .collect();

    let mut h = vec![0.0f64; ndof * ndof];
    for dof in 0..ndof {
        let gf = &grads[2 * dof];
        let gb = &grads[2 * dof + 1];
        for j in 0..ndof {
            h[dof * ndof + j] = (gf[j] - gb[j]) / (2.0 * step_size);
        }
    }

    // Symmetrize H_ij ← (H_ij + H_ji) / 2
    for i in 0..ndof {
        for j in i + 1..ndof {
            let avg = (h[i * ndof + j] + h[j * ndof + i]) / 2.0;
            h[i * ndof + j] = avg;
            h[j * ndof + i] = avg;
        }
    }

    h
}

/// Clone `molecule` with atom `idx` displaced by `delta` along `axis` (0=x, 1=y, 2=z).
pub fn displaced(molecule: &Molecule, idx: usize, axis: usize, delta: f64) -> Molecule {
    let mut atoms = molecule.atoms.clone();
    atoms[idx].position[axis] += delta;
    Molecule::new(atoms, molecule.charge, molecule.multiplicity)
}