Function initialize_inverse_hessian 

pub fn initialize_inverse_hessian(n: usize) -> DMatrix<f64>

Updates the Hessian matrix using the PSB (Powell-Symmetric-Broyden) formula.

The PSB formula is a rank-2 update that modifies the Hessian approximation based on the difference in gradients (yk) and the step taken (sk):

H_new = H + (yk - H*sk) * sk^T + sk * (yk - H*sk)^T
        - [(yk - H*sk)^T * sk] * (sk * sk^T) / (sk^T * sk)

This update preserves symmetry and positive definiteness under certain conditions. The PSB update is more stable than BFGS for poorly conditioned problems.

Arguments

• hessian - Current Hessian approximation
• sk - Step vector (x_new - x_old)
• yk - Gradient difference (g_new - g_old)

Returns

Returns the updated Hessian matrix as a DMatrix<f64>.

Examples

use omecp::optimizer::update_hessian_psb;
use nalgebra::{DMatrix, DVector};

let h_old = DMatrix::identity(3, 3);
let sk = DVector::from_vec(vec![0.1, 0.2, 0.3]);
let yk = DVector::from_vec(vec![0.05, 0.1, 0.15]);

// let h_new = update_hessian_psb(&h_old, &sk, &yk);

Initializes the inverse Hessian matrix for BFGS optimization.

Following the Fortran MECP implementation (adapted for Angstrom), the inverse Hessian is initialized as a diagonal matrix with value 0.7 Ų/Ha. This corresponds to a Hessian of approximately 1.4 Ha/Ų, which provides reasonable initial step sizes.

Arguments

• n - Dimension of the matrix (3 × number of atoms)

Returns

Returns an n×n diagonal matrix with 0.7 on the diagonal.

Units

The inverse Hessian is in Ų/Ha (Angstrom squared per Hartree). This matches the Angstrom-based unit system used throughout OpenMECP:

• Newton step: H⁻¹ (Ų/Ha) × g (Ha/A) = step (A)
• The 0.7 value provides reasonable initial step sizes for molecular systems

Requirements

Validates: Requirements 3.1, 8.2