Function bfgs_step 

pub fn bfgs_step(
    x0: &DVector<f64>,
    g0: &DVector<f64>,
    inv_hessian: &DMatrix<f64>,
    config: &Config,
    _adaptive_scale: f64,
) -> DVector<f64>

Performs a BFGS optimization step.

BFGS (Broyden-Fletcher-Goldfarb-Shanno) is a quasi-Newton optimization method that approximates the inverse Hessian using gradient information. It provides good convergence for the first few iterations while building curvature information.

The BFGS step direction is computed by solving:

d = -H^(-1) * g

where H is the Hessian approximation and g is the gradient. The step size is automatically limited by config.max_step_size to prevent overshooting.

Arguments

• x0 - Current geometry coordinates
• g0 - Current MECP gradient
• hessian - Current Hessian approximation matrix
• config - Configuration with step size limits and other parameters

Returns

Returns the new geometry coordinates after the BFGS step as a DVector<f64>.

Examples

use omecp::optimizer::bfgs_step;
use nalgebra::DVector;

let x0 = DVector::from_vec(vec![0.0, 0.0, 0.0]);
let g0 = DVector::from_vec(vec![0.1, 0.2, 0.3]);
let hessian = DMatrix::identity(3, 3);

// let x_new = bfgs_step(&x0, &g0, &hessian, &config, 1.0);

pub fn bfgs_step( x0: &DVector, g0: &DVector, hessian: &DMatrix, config: &Config, _adaptive_scale: f64, // Parameter kept for compatibility but not used for BFGS ) -> DVector { // Exact propagationBFGS implementation: // 1. dk = -H^-1 * g (Newton direction) // 2. if ||dk|| > 0.1: dk = dk * 0.1 / ||dk|| (cap direction to 0.1 Angstrom) // 3. XNew = X0 + rho * dk (rho=15 for MECP) // 4. MaxStep: if ||XNew - X0|| > MAX_STEP_SIZE: scale to MAX_STEP_SIZE

// Step 1: Compute Newton direction dk = -H^-1 * g
let neg_g = -g0;
let mut dk = hessian.clone().lu().solve(&neg_g).unwrap_or_else(|| {
    // Fallback to steepest descent when Hessian is singular
    println!("BFGS Step: Hessian is singular, falling back to steepest descent");
    -g0 / (g0.norm() + 1e-14)
});

// Step 2: Cap dk to 0.1 Angstrom 
// Convert to Bohr since internal coordinates are in Bohr
let dk_cap = 0.1 * ANGSTROM_TO_BOHR; // 0.1 Angstrom in Bohr
let dk_norm = dk.norm();
if dk_norm > dk_cap {
    println!(
        "BFGS: dk norm {:.6} > {:.6}, capping direction",
        dk_norm, dk_cap
    );
    dk *= dk_cap / dk_norm;
}

// Step 3: Apply rho multiplier (rho=15 for MECP optimization)
// This aggressive multiplier helps escape shallow regions quickly
// Note: dk is in Bohr (same as coordinates), so no unit conversion needed
let rho = config.bfgs_rho;
let x_new = x0 + &dk * rho;

// Step 4: MaxStep - limit total step to max_step_size
let step = &x_new - x0;
let step_norm = step.norm();

// Debug: print step details
let step_angstrom = step_norm * crate::config::BOHR_TO_ANGSTROM;
println!(
    "BFGS: dk_norm={:.6}, dk_capped={:.6}, rho={}, raw_step={:.6} bohr ({:.6} Ang)",
    dk_norm, dk.norm(), rho, step_norm, step_angstrom
);

if step_norm > config.max_step_size {
    let scale = config.max_step_size / step_norm;
    let final_step_angstrom = config.max_step_size * crate::config::BOHR_TO_ANGSTROM;
    println!(
        "BFGS step: {:.6} -> {:.6} bohr ({:.6} Ang) (MaxStep applied)",
        step_norm, config.max_step_size, final_step_angstrom
    );
    x0 + &step * scale
} else {
    println!("BFGS step: {:.6} bohr ({:.6} Ang) (within max_step_size)", step_norm, step_angstrom);
    x_new
}

} Performs a BFGS optimization step following the Fortran MECP implementation.

This implementation is adapted from the original Harvey Fortran code (eas) but operates in Angstrom-based units:

• Uses inverse Hessian (Ų/Ha) for Newton step
• Works in Angstrom for the Newton step computation
• Two-stage step limiting: total norm, then max component (in A)

Algorithm (from Fortran UpdateX subroutine, adapted for Angstrom)

1. First step: ChgeX = -0.7 * G (steepest descent with H_inv diagonal = 0.7)
2. Later steps: ChgeX = -H_inv * G (Newton step with BFGS-updated inverse Hessian)
3. Limit step: if ||ChgeX|| > 0.1*N A, scale down
4. Limit components: if max(|ChgeX_i|) > 0.1 A, scale down
5. Add Angstrom step to Angstrom coordinates

Units

• Input coordinates (x0): Angstrom
• Input gradient (g0): Ha/A (converted from native QM output)
• Inverse Hessian: Ų/Ha (initialized to 0.7 on diagonal)
• Newton step: A (H⁻¹ × g = Ų/Ha × Ha/A = A)
• Output coordinates: Angstrom

Requirements

Validates: Requirements 4.1, 4.2, 4.3, 4.4, 4.5