Block-Jacobi preconditioner for EwaldPeriodicDipoleOperator's own PCG solve.
USAGE: Eigen::ConjugateGradient owns its own preconditioner as a plain (default-constructed) member (Eigen::IterativeSolverBase<>:: m_preconditioner), and cg.compute(op) calls that member's own compute(op.matrix()) – it never passes a caller-supplied preconditioner instance through. Confirmed directly (a real, not hypothetical, failure caught by an actual compiled-and-run test): calling EwaldBlockJacobiPreconditioner(registry, ids, alpha, thole_a) and then handing that instance to a template parameter alone is NOT enough – cg.compute(op) still default-constructs its own (uninitialized) instance over it. The construct that actually works, verified end to end via Eigen::ConjugateGradient with a genuine block-diagonal test operator (0 iterations to machine precision, since the preconditioner exactly matched the operator there):
Eigen::ConjugateGradient<Op, Lower|Upper, EwaldBlockJacobiPreconditioner> cg; cg.compute(op); cg.preconditioner() = EwaldBlockJacobiPreconditioner(registry, ids, alpha_ewald, thole_a); x = cg.solveWithGuess(b, x0);
i.e. assign the real, already-built preconditioner via IterativeSolverBase<>::preconditioner()'s own public, non-const accessor, AFTER compute(), not by relying on compute() to build it.
This exists because Eigen::DiagonalPreconditioner<double> – the default used with EwaldPeriodicDipoleOperator – reads only operator()(i,j)'s own same-site diagonal (getPInv()), with zero knowledge of the intramolecular (Thole-damped) coupling EwaldPeriodicDipoleOperator:: AddIntraSegmentCoupling actually adds for multi-site segments. That gap was found, this session, to matter a great deal in practice: enabling Thole damping on the intramolecular term (the physically correct, legacy-matching behavior – confirmed directly against legacy's own FU12_ERFC_At_By/UpdateAllBls) made a real, large ~5000-site PCG solve diverge outright (residual growing past its own starting value, not merely converging slowly) with the bare-diagonal preconditioner, while disabling that same damping (reverting to the undamped erfc-only tensor) let the identical system converge cleanly in 51 iterations. Since the physics itself is confirmed correct (matching legacy), the fix pursued here is a better preconditioner, not reverting the physics.
A direct numerical check (methane geometry, realistic atomic polarizabilities, thole_a=0.6) found the intramolecular block ALONE, isolated from the rest of the system, IS positive-definite when Thole-damped (min eigenvalue ~0.05, comfortably positive) and is NOT when undamped (min eigenvalue ~-0.12) – the opposite of what a first, trace-based argument suggested (a real correction made mid-session: a negative trace does not by itself imply a negative eigenvalue exists, and here it doesn't). This means the observed divergence is NOT because the local, per-segment intramolecular coupling is itself indefinite – whatever is happening lives either in how that locally-SPD block combines with the periodic (intermolecular, shape, self-field) terms across the full ~5000-site system, or is purely a preconditioner- tracking problem (the true global operator staying SPD throughout, but DiagonalPreconditioner's blindness to the intramolecular term causing CG's implicit preconditioned-residual bookkeeping to behave as though it doesn't). This class was not built having resolved which of those is true – only having confirmed the local block itself is a valid (SPD) building block either way, and that a preconditioner genuinely reflecting the intramolecular coupling should track the true operator far better than the bare diagonal does regardless of which explanation is correct. Whether this actually resolves the real ~5000-site divergence, as opposed to just being well-motivated, has not yet been tested against that real system – only the mechanics above have been verified, with a synthetic operator, not the real EwaldPeriodicDipole Operator/real molecular system.
What this computes: for every segment in ids with 2+ sites, the full dense 3*N_site x 3*N_site local block – getPInv() on each site's own diagonal, PLUS the same Thole-damped intramolecular coupling AddIntraSegmentCoupling itself adds (identical B-function/ComputeThole calls, same sign convention) – factorized ONCE via LDLT (chosen over LLT/Cholesky specifically because the local block was found to be only marginally SPD at higher thole_a, e.g. min eigenvalue ~0.0016 at thole_a=1.0 in the same methane test above – LDLT degrades gracefully on a near-singular or, for some future geometry/polarizability combination, genuinely indefinite block, where LLT would simply fail). Single-site segments fall back to the plain diagonal (getPInv()), identical to what DiagonalPreconditioner itself would do for them.
solve(b) applies each segment's own precomputed factorization to that segment's own slice of b independently – a genuine block-Jacobi scheme: exact within each segment's own local block, still ignoring inter-segment (periodic, real+reciprocal-space) coupling entirely, the same simplification DiagonalPreconditioner itself already makes at the single-site level. The intermolecular coupling is left to PCG's own outer iteration to resolve, same as before – only the intramolecular piece moves from "completely ignored" to "handled exactly."
Definition at line 110 of file ewaldblockjacobipreconditioner.h.
