Description:

This task represents the final core implementation phase before exposing the public API. It focuses on creating the QEqSolver engine, which orchestrates the entire charge equilibration process. This solver will integrate the atomic parameters from the params module and the shielded Coulomb potential from the math module to iteratively solve for equilibrium charges.

The implementation will reside in a new top-level solver module and will use the faer crate for high-performance, numerically stable linear algebra. The core logic will be a Self-Consistent Field (SCF) loop to handle the charge-dependent hardness of hydrogen, making the solver applicable to all elements described in the Rappe & Goddard paper.

Tasks:

• Phase 1: Module Scaffolding and Dependencies

  • Add the faer crate (version = "0.17") to Cargo.toml.
  • Create a new directory src/solver/ with mod.rs.
  • Create src/solver/options.rs and define the SolverOptions struct with tolerance, max_iterations, and lambda_scale, including a Default implementation.
  • Create src/solver/solver.rs and define the QEqSolver struct. It should hold a reference to Parameters and an instance of SolverOptions.
  • Implement QEqSolver::new() and a builder-style with_options() method.

• Phase 2: Implement the System Builder (build_system)

  • In solver.rs, create the private helper function build_system.
  • Signature: fn build_system<A: AtomView>(&self, atoms: &[A], element_data: &[&'p ElementData], current_charges: ColRef<f64>) -> Result<(Mat<f64>, Col<f64>), CheqError>.
  • Logic:
    • Initialize an (N+1) x (N+1) matrix A and an (N+1) vector b using faer.
    • Loop i from 0 to N-1:
      • Diagonal A[i, i] (Hardness J_ii):
        • Check if the atom is Hydrogen (n=1 or atomic_number=1).
        • If yes, apply the charge-dependent hardness formula: J_ii = J_ii⁰ * (1 + Q_i * H_CHARGE_FACTOR).
        • If no, use the hardness J_ii⁰ directly from element_data.
      • Off-diagonal A[i, j] (Coulomb J_ij):
        • Loop j from i+1 to N-1.
        • Calculate inter-atomic distance R_ij in Angstroms.
        • Perform all necessary unit conversions (Å → Bohr) for inputs to the math module.
        • Call math::shielding::gaussian_coulomb_integral.
        • Convert the result (Hartree → eV).
        • Fill the symmetric elements A[i, j] and A[j, i].
      • Right-hand side b[i] (Electronegativity χ_i⁰):
        • Set b[i] = -element_data[i].electronegativity.
    • Constraint Equations:
      • Fill the last column of A (rows 0..N) with -1.0.
      • Fill the last row of A (cols 0..N) with 1.0. A[N, N] remains 0.0.
      • Set the last element of the vector b[N] = total_charge.
  • Return Ok((A, b)).

• Phase 3: Implement the SCF Solver Logic (solve)

  • Implement the public method pub fn solve<A: AtomView>(&self, atoms: &[A], total_charge: f64) -> Result<CalculationResult, CheqError>.
  • Initial Setup:
    • Handle the edge case of an empty atoms slice, returning CheqError::NoAtoms.
    • Pre-fetch all necessary ElementData for the atoms, handling CheqError::ParameterNotFound.
    • Initialize charges vector to all zeros.
    • Determine if the system contains hydrogen to decide if an SCF loop is needed.
  • SCF Loop:
    • Loop up to self.options.max_iterations.
    • Call build_system to get the current matrix A and vector b.
    • Solve Ax = b using a.lu().solve(&b). Handle potential None result by returning CheqError::LinalgError.
    • Extract new charges Q_new from the first N elements of the solution vector.
    • Calculate max_charge_delta = max(|Q_new - Q_old|).
    • Update charges = new_charges.
    • Check for Convergence: If !has_hydrogen (exit after one iteration) or max_charge_delta < self.options.tolerance, break the loop and proceed to success.
  • Finalization:
    • If the loop converged, extract the equilibrated potential χ̄ = -solution[N], create CalculationResult, and return Ok(...).
    • If the loop terminated due to max iterations, return Err(CheqError::NotConverged).

• Phase 4: Integration Testing

  • Create a tests submodule within src/solver/solver.rs.
  • Test a linear system (no H): Use a molecule like CO or N₂. Assert that charges are computed in 1 iteration and have the correct signs and rough magnitudes.
  • Test a non-linear system (with H):
    • Use LiH. Check that it takes multiple iterations. Assert the final charge on H is negative and close to the paper's value (~ -0.7).
    • Use H₂O. Assert the final charge on H is positive and close to the paper's QEq value (~ 0.35).
  • Test symmetry: For H₂, assert charges are exactly zero. For a symmetric molecule like CH₄, assert all H charges are identical.
  • Test error handling:
    • Write a test that calls solve with an element not in the default parameters (e.g., Uranium, if not present) and assert it returns CheqError::ParameterNotFound.
    • Write a test with an empty slice of atoms and assert it returns CheqError::NoAtoms.