Function bfgs_step_direct 

pub fn bfgs_step_direct(
    x0: &DVector<f64>,
    g0: &DVector<f64>,
    hessian: &DMatrix<f64>,
    config: &Config,
) -> DVector<f64>

Performs a BFGS step using a direct Hessian (matching exactly).

This is a direct port of propagationBFGS() + MaxStep().

Algorithm )

1. dk = solve(Bk, -Gk)              # Newton direction via LU decomposition
2. if ||dk|| > CAP: dk *= CAP/||dk||  # Cap direction magnitude
3. step = rho * dk                    # Amplify small Newton steps
4. if ||step|| > MAX: step *= MAX/||step||  # Final MaxStep cap
5. XNew = X0 + step

Arguments

• x0 - Current geometry coordinates in A
• g0 - Current MECP gradient (mixed units: Ha + Ha/A)
• hessian - Current direct Hessian matrix (Ha/A²)
• config - Configuration with step size limits

Returns

Returns the new geometry coordinates in A.

Units

• dk = B⁻¹ × g: (Ha/A²)⁻¹ × (Ha/A) = A (step in Angstrom)
• Cap: 0.1 (NOT in Bohr, operates in mixed-unit Newton space)
• rho: 15.0 (amplification factor)
• MaxStep: config.max_step_size in A (default: 0.1 A)