Function build_gediis_b_matrix 

fn build_gediis_b_matrix(
    grads: &VecDeque<DVector<f64>>,
    f_vecs: &VecDeque<DVector<f64>>,
    hessians: &VecDeque<DMatrix<f64>>,
    energies: &VecDeque<f64>,
) -> DMatrix<f64>

Computes enhanced error vectors for GEDIIS optimization.

GEDIIS error vectors incorporate both gradient and energy information to provide better convergence for MECP optimization. The energy contribution helps emphasize geometries that are closer to the target energy difference.

Arguments

• grads - History of gradient vectors from previous iterations
• energies - History of energy differences (E1 - E2) from previous iterations

Returns

Returns a vector of enhanced error vectors that include energy weighting. Each error vector combines gradient information with energy deviation.

Algorithm

For each iteration i:

error[i] = g[i] + λ * (E[i] - E_avg) * g[i]

where:

• g[i] is the gradient at iteration i
• E[i] is the energy difference at iteration i
• E_avg is the average energy difference over all iterations
• λ = 0.05 is a FIXED small constant (typically 0.01-0.1)

Important: Fixed Lambda

The lambda parameter MUST be fixed and small (0.01-0.1), NOT adaptive. Using adaptive scaling like λ = 0.1/|g| causes catastrophic instability near convergence because:

• When |g| → 0, λ → ∞
• Tiny energy noise (10⁻⁸) gets amplified to 10⁻¹ in error vector
• Destroys convergence

Reference: Truhlar et al., J. Chem. Theory Comput. 2006, 2, 835-839 explicitly warns against adaptive scaling.

Builds the B matrix for standard GEDIIS optimization.

Uses the formula from Li, Frisch, and Truhlar (J. Chem. Theory Comput. 2006, 2, 835-839):

B[i,j] = -(g_i - g_j) · (x_i - x_j)

This metric captures the curvature of the energy surface without explicit Hessian.

Unit Analysis

• g_i - g_j: Gradient difference in Ha/A
• x_i - x_j: Geometry difference in Angstrom
• B[i,j]: Mixed units (Ha/A × A = Ha)

The mixed units are acceptable because the B-matrix is only used to solve for dimensionless interpolation coefficients. The DIIS constraint (Σc_i = 1) ensures the coefficients are scale-invariant.

Arguments

• grads - History of gradient vectors in Ha/A
• geoms - History of geometry vectors in Angstrom

Returns

Returns the (n+1) × (n+1) B matrix for DIIS coefficient determination. Builds a stable GEDIIS B-matrix using GDIIS-style error vectors with energy coupling.

B[i,j] = e_i·e_j + α·E_i·E_j

where:

• e_i = H̄⁻¹ · (g_i + f_i): Newton-step error vectors (same as GDIIS)
• E_i = energy gap at point i (MECP condition)
• α = mean(|e·e|) / mean(|E·E|): dynamically balanced coupling

Compared to the old formulation -(g_i-g_j)·(x_i-x_j) which was ill-conditioned (all entries tiny and nearly identical), this uses the well-conditioned GDIIS error vectors with a small energy bias.