Description:

This task covers the implementation of the mathematical core for charge interaction in the cheq library, located in the math::shielding module. The chosen model is a numerically robust, Gaussian-screened Coulomb potential, which faithfully approximates the behavior of Slater-Type Orbitals (STOs) while avoiding the complexity of direct STO integration. This is effectively an STO-1G (Slater-Type Orbital approximated by 1 Gaussian) approach.

The primary function, gaussian_coulomb_integral, takes fundamental atomic parameters (n, r_cov_bohr, lambda) and the inter-atomic distance in Bohr as input. It returns the shielded electrostatic interaction energy (J_AB) in Hartree atomic units. This function is the cornerstone for building the off-diagonal elements of the QEq matrix.

The implementation is self-contained within the math::shielding module, deriving all necessary intermediate parameters (like Slater ζ and equivalent Gaussian α exponents) internally. It correctly handles the R → 0 limit case to ensure physical behavior and numerical stability.

Tasks:

• Phase 1: Module and Function Scaffolding

  • Create the file src/math/shielding.rs.
  • Define the public function signature: pub fn gaussian_coulomb_integral(distance_bohr: f64, n1: u8, r_cov1_bohr: f64, n2: u8, r_cov2_bohr: f64, lambda: f64) -> f64.
  • Document the function, specifying that all inputs are in atomic units (Bohr) and the output is in Hartree.
  • Define necessary private helper functions: slater_exponent, equivalent_gaussian_exponent, get_c_n_fit_coeff, and screened_potential.
  • Define required physical constants in src/math/constants.rs (BOHR_TO_ANGSTROM, HARTREE_TO_EV).

• Phase 2: Implement the Gaussian Screening Algorithm

  • Inside gaussian_coulomb_integral:
    • Step 1: Calculate Slater Exponents (ζ)
      • Call slater_exponent for both atoms, using their n, r_cov_bohr, and the global lambda. The formula ζ = λ * (2n + 1) / (2 * R_bohr) is implemented.
    • Step 2: Calculate Equivalent Gaussian Exponents (α)
      • Call equivalent_gaussian_exponent for both atoms, using their n and calculated ζ. The formula α = c(n) * ζ² / n is implemented, with c(n) values from get_c_n_fit_coeff.
    • Step 3: Calculate the Screened Potential
      • Call the screened_potential helper function.
      • Internally, it calculates the screening parameter β = sqrt(2 * α_1 * α_2 / (α_1 + α_2)).
      • It then calculates the shielded interaction in Hartree atomic units: J_hartree = erf(β * R_bohr) / R_bohr. This includes the special handling for the R → 0 limit using the Taylor expansion of erf.
    • Step 4: Return Result
      • Return the calculated J_hartree. The conversion to eV is handled by the caller (the solver) to keep the math module in pure atomic units.

• Phase 3: Rigorous Unit Testing

  • In src/math/shielding.rs, create a tests submodule.
  • Test the R → ∞ asymptotic behavior:
    • A test test_integral_at_large_distance_asymptote asserts that for a large distance (R = 1000.0 Bohr), the result is very close to the point-charge Coulomb law 1.0 / R_bohr.
  • Test specific known values:
    • A test test_integral_with_known_values verifies calculated J_AB values for H₂ and a C-O pair against pre-computed benchmarks to ensure correctness.
  • Test numerical stability near R → 0:
    • A test test_integral_at_zero_distance_limit asserts that the function returns a finite, positive value at R = 0 and that the value is consistent with the limit as R approaches zero, confirming the erf expansion is working correctly.
  • Test symmetry:
    • A test test_integral_symmetry_property asserts that J(A, B) is identical to J(B, A).
  • Test chemical and physical trends:
    • A test test_integral_chemical_trends confirms that the interaction decreases with distance and that interactions between more diffuse orbitals (e.g., Na-Na) are weaker than between more compact orbitals (e.g., H-H) at the same separation.
    • A test test_integral_scaling_with_lambda confirms the physical expectation that a smaller lambda leads to stronger screening (smaller J_AB).