chemx-props 0.1.0

Molecular properties and frequencies for chemx: dipole, populations, bond orders, numerical Hessian, harmonic frequencies, RRHO thermochemistry.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

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

use chemx_core::Molecule;

/// 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.
pub fn numerical_hessian<F>(molecule: &Molecule, step_size: f64, gradient_fn: F) -> Vec<f64>
where
    F: Fn(&Molecule) -> Vec<f64>,
{
    let natom = molecule.len();
    let ndof = 3 * natom;
    let mut h = vec![0.0f64; ndof * ndof];

    for atom in 0..natom {
        for axis in 0..3 {
            let dof = atom * 3 + axis;
            let mol_fwd = displaced(molecule, atom, axis, step_size);
            let mol_bwd = displaced(molecule, atom, axis, -step_size);
            let gf = gradient_fn(&mol_fwd);
            let gb = gradient_fn(&mol_bwd);
            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)
}