Function compute_mecp_gradient 

pub fn compute_mecp_gradient(
    state_a: &State,
    state_b: &State,
    fixed_atoms: &[usize],
) -> MecpGradient

Computes the MECP effective gradient for optimization.

This function implements the Harvey et al. algorithm for MECP optimization by computing the effective gradient that drives the system toward the minimum energy crossing point. The gradient has two components:

1. f-vector: Drives the energy difference (E1 - E2) to zero
2. g-vector: Minimizes the energy perpendicular to the gradient difference

The effective gradient is computed as:

G = (E1 - E2) * x_norm + [f1 - (x_norm · f1) * x_norm]
    \_____f-vector____/   \________g-vector________/

where x_norm = (f1 - f2) / |f1 - f2| is the normalized gradient difference.

Unit Analysis

This implementation matches the reference (adapted for Angstrom units):

• Input forces (state_a.forces, state_b.forces): Ha/A (converted from native QM output)
• f1, f2 (negated forces = gradients): Ha/A
• x_norm (normalized gradient difference): dimensionless (unit vector)
• f-vector = (E1 - E2) × x_norm: Hartree (energy × dimensionless)
• g-vector = f1 - (x_norm · f1) × x_norm: Ha/A
• Combined gradient: Mixed units (Ha + Ha/A)

The mixed units are intentional in the Harvey algorithm:

• The f-vector acts as a penalty term driving energy difference to zero
• The g-vector minimizes energy perpendicular to the crossing seam
• Both components contribute appropriately to the optimization direction

When used with the BFGS optimizer, the inverse Hessian (Ų/Ha) handles the unit conversion to produce steps in Angstrom.

Arguments

• state_a - Electronic state 1 (energy in Ha, forces in Ha/A, geometry)
• state_b - Electronic state 2 (energy in Ha, forces in Ha/A, geometry)
• fixed_atoms - List of atom indices to fix during optimization (0-based)

Returns

Returns the MECP effective gradient as a DVector<f64> with length 3 × num_atoms. The gradient has mixed units (f-vector in Ha, g-vector in Ha/A) matching the reference implementation.

Requirements

Validates: Requirements 6.1, 6.2, 6.3

Examples

use omecp::geometry::{Geometry, State};
use omecp::optimizer::compute_mecp_gradient;

// let gradient = compute_mecp_gradient(&state_a, &state_b, &[]);
// assert_eq!(gradient.len(), state_a.geometry.num_atoms * 3);