chemx-props 0.3.0

//! Single-Point Hessian (SPH) vibrational analysis.
//!
//! Reference: S. Spicher and S. Grimme, "Single-Point Hessian Calculations for
//! Improved Vibrational Frequencies and Rigid-Rotor-Harmonic-Oscillator
//! Thermostatistics", *J. Chem. Theory Comput.* 2021, **17**, 1701–1714
//! (DOI 10.1021/acs.jctc.0c01306).
//!
//! ## What this implements
//!
//! Ordinary harmonic frequency analysis assumes the Hessian is evaluated at a
//! stationary point (gradient ≈ 0). The SPH idea is to extract meaningful
//! RRHO frequencies from a Hessian computed at a **non-stationary** geometry
//! (the structure is taken as-is, with no preceding optimization), by removing
//! the contribution of the non-zero gradient.
//!
//! chemx implements the **gradient-projection** form of this treatment: in
//! addition to the six (five for a linear molecule) Eckart translation/rotation
//! vectors, the normalized mass-weighted gradient direction is projected out of
//! the mass-weighted Hessian before diagonalization. This is the reaction-path
//! projection of Miller, Handy & Adams (*J. Chem. Phys.* 1980, **72**, 99) — the
//! gradient direction is exactly the steepest-descent (reaction-path) mode at a
//! non-stationary geometry, and removing it suppresses the spurious soft/
//! imaginary mode it would otherwise mix into the spectrum. In mass-weighted
//! coordinates `q_i = √m_i x_i`, the gradient transforms as
//! `g^mw_i = g_i / √m_i`.
//!
//! **Key property (and the validity test):** at a true minimum the gradient is
//! ≈ 0, its mass-weighted vector has ≈ 0 norm and is dropped by the
//! Gram–Schmidt step, so the SPH projector reduces *exactly* to the ordinary
//! Eckart projector and the SPH frequencies reproduce the ordinary harmonic
//! frequencies.
//!
//! ## Deviation from xtb `--bhess`
//!
//! This is **not** identical to xtb's `--bhess`, which instead adds a
//! root-mean-square-deviation (RMSD) Gaussian biasing potential (the xtb
//! metadynamics bias, with method-specific `kpush`/`alpha`) centred on the
//! input geometry and computes the Hessian of that biased surface. The
//! gradient-projection form here is a documented, parameter-free alternative
//! with the same intent and the same minimum-reproduction property; it does not
//! introduce the bias-strength parameters and so makes no claim to reproduce
//! xtb `--bhess` numbers. See `PROVENANCE.md`. SPH frequencies are approximate
//! and the driver labels them as such.

use chemx_core::Molecule;

use crate::frequencies::{FrequencyResult, harmonic_frequencies_projected};

/// Threshold on the mass-weighted gradient 2-norm below which the geometry is
/// treated as stationary and the gradient projection is a no-op — the SPH
/// result then equals the ordinary harmonic analysis (all 3N−6/3N−5 modes
/// retained). The gradient direction, once normalized, removes a full mode, so
/// the reduction to ordinary frequencies must be a clean on/off decision rather
/// than an O(|g|) effect. The value sits safely between a tight-optimization
/// residual (mass-weighted norm ~1e-5 for chemx's `max_force = 3e-6` Eh/bohr
/// tight criterion) and any chemically displaced gradient (~1e-2 and up).
pub const SPH_GRADIENT_NORM_FLOOR: f64 = 1.0e-3;

/// Compute SPH harmonic frequencies from a Cartesian Hessian and the Cartesian
/// gradient at the **same (non-stationary) geometry**.
///
/// * `hessian` — 3N×3N Cartesian Hessian (row-major), hartree/bohr².
/// * `gradient_cart` — flat Cartesian gradient `[gx₁,gy₁,gz₁,…]` (length 3N),
///   hartree/bohr, at the geometry of `molecule`.
///
/// Builds the mass-weighted gradient direction and projects it out alongside
/// translations/rotations. If the gradient norm is below
/// [`SPH_GRADIENT_NORM_FLOOR`], no gradient projection is applied and the result
/// equals [`crate::frequencies::harmonic_frequencies`].
pub fn sph_frequencies(
    molecule: &Molecule,
    hessian: &[f64],
    gradient_cart: &[f64],
) -> FrequencyResult {
    let natom = molecule.len();
    let ndof = 3 * natom;
    assert_eq!(gradient_cart.len(), ndof, "gradient must be length 3N");

    // Mass-weight the gradient: g^mw_i = g_i / √m_i.
    let masses: Vec<f64> = molecule.atoms.iter().map(|a| a.element.mass()).collect();
    let mut g_mw = vec![0.0f64; ndof];
    for i in 0..natom {
        let sm = masses[i].sqrt();
        for k in 0..3 {
            g_mw[3 * i + k] = gradient_cart[3 * i + k] / sm;
        }
    }
    let norm: f64 = g_mw.iter().map(|x| x * x).sum::<f64>().sqrt();

    if norm < SPH_GRADIENT_NORM_FLOOR {
        // Stationary (or numerically so): plain Eckart projection.
        harmonic_frequencies_projected(molecule, hessian, &[])
    } else {
        // Normalize and project out the gradient direction in addition.
        for x in &mut g_mw {
            *x /= norm;
        }
        harmonic_frequencies_projected(molecule, hessian, &[g_mw])
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::frequencies::harmonic_frequencies;
    use crate::hessian::numerical_hessian;
    use chemx_core::{Atom, Element};

    /// A toy quadratic surface E(x) = ½ (x−x₀)ᵀ K (x−x₀) + gᵀ(x−x₀) lets us test
    /// the projector analytically without an SCF: its gradient is K(x−x₀)+g and
    /// its Hessian is the constant K. We build a diatomic-like 2-atom system.
    fn toy() -> (Molecule, Vec<f64>) {
        // Two atoms on the z-axis; a simple spring along z plus stiff x/y so the
        // Hessian is well-conditioned. Masses from real elements (H, H).
        let mol = Molecule::new(
            vec![
                Atom::new(Element::from_z(1).unwrap(), [0.0, 0.0, -0.7]),
                Atom::new(Element::from_z(1).unwrap(), [0.0, 0.0, 0.7]),
            ],
            0,
            1,
        );
        // A 6×6 Hessian: stretch along z with force constant 0.5, plus weak
        // restoring terms so trans/rot project cleanly. Use a physically shaped
        // Hessian: k for the symmetric stretch coordinate.
        let k = 0.5;
        let mut h = vec![0.0f64; 36];
        // Couple the two z coordinates as a spring: H_zz = [[k,-k],[-k,k]].
        let zz = [(2usize, 2usize, k), (2, 5, -k), (5, 2, -k), (5, 5, k)];
        for (i, j, v) in zz {
            h[i * 6 + j] = v;
        }
        (mol, h)
    }

    #[test]
    fn reduces_to_ordinary_at_zero_gradient() {
        let (mol, h) = toy();
        let zero_grad = vec![0.0f64; 6];
        let sph = sph_frequencies(&mol, &h, &zero_grad);
        let ord = harmonic_frequencies(&mol, &h);
        assert_eq!(sph.frequencies_cm1.len(), ord.frequencies_cm1.len());
        for (a, b) in sph.frequencies_cm1.iter().zip(&ord.frequencies_cm1) {
            assert!((a - b).abs() < 1e-9, "SPH {a} vs ordinary {b}");
        }
    }

    #[test]
    fn tiny_gradient_below_floor_is_noop() {
        let (mol, h) = toy();
        let grad = vec![1e-9; 6]; // norm below floor
        let sph = sph_frequencies(&mol, &h, &grad);
        let ord = harmonic_frequencies(&mol, &h);
        for (a, b) in sph.frequencies_cm1.iter().zip(&ord.frequencies_cm1) {
            assert!((a - b).abs() < 1e-9);
        }
    }

    #[test]
    fn gradient_projection_removes_a_mode() {
        // With a non-zero gradient, one additional direction is projected out,
        // so the SPH spectrum has one more ~zero eigenvalue than the ordinary
        // one (the gradient mode is removed from the vibrational block).
        let (mol, h) = toy();
        let grad = vec![0.0, 0.0, 0.3, 0.0, 0.0, -0.3]; // along the stretch
        let sph = sph_frequencies(&mol, &h, &grad);
        // The stretch mode (the only genuine vibration of this toy) is the
        // gradient direction, so projecting it out leaves no real vibration:
        // all modes are ~0 (trans/rot/gradient) and none are imaginary.
        assert_eq!(sph.n_imaginary, 0);
        let n_nonzero = sph
            .frequencies_cm1
            .iter()
            .filter(|f| f.abs() > 10.0)
            .count();
        let ord = harmonic_frequencies(&mol, &h);
        let n_nonzero_ord = ord
            .frequencies_cm1
            .iter()
            .filter(|f| f.abs() > 10.0)
            .count();
        assert_eq!(n_nonzero_ord, 1, "toy has 1 vibration");
        assert_eq!(n_nonzero, 0, "gradient mode projected out");
    }

    #[test]
    fn sph_water_displaced_no_spurious_imaginary() {
        // A real (small) check that the machinery is NaN-free and that a
        // displaced water Hessian built by finite differences of a toy harmonic
        // gradient yields no spurious imaginary modes after SPH projection.
        // (The end-to-end SCF version lives in chemx integration tests.)
        let mol = Molecule::new(
            vec![
                Atom::new(Element::from_z(8).unwrap(), [0.0, 0.0, 0.20]),
                Atom::new(Element::from_z(1).unwrap(), [0.0, 1.43, -0.90]),
                Atom::new(Element::from_z(1).unwrap(), [0.0, -1.43, -0.90]),
            ],
            0,
            1,
        );
        // Harmonic model around a reference; gradient is non-zero at `mol`.
        let ref_pos: Vec<[f64; 3]> =
            vec![[0.0, 0.0, 0.12], [0.0, 1.43, -0.95], [0.0, -1.43, -0.95]];
        let k = 0.4;
        let grad_fn = |m: &Molecule| -> Vec<f64> {
            let mut g = vec![0.0f64; 9];
            for i in 0..3 {
                for c in 0..3 {
                    g[3 * i + c] = k * (m.atoms[i].position[c] - ref_pos[i][c]);
                }
            }
            g
        };
        let hess = numerical_hessian(&mol, 0.005, |m| grad_fn(m));
        let grad_here = grad_fn(&mol);
        let sph = sph_frequencies(&mol, &hess, &grad_here);
        assert!(sph.frequencies_cm1.iter().all(|f| f.is_finite()));
        assert_eq!(
            sph.n_imaginary, 0,
            "spurious imaginary: {:?}",
            sph.frequencies_cm1
        );
    }
}