Bug: the “cluster” sub-calculation lands on the wrong SCF solution near the Na⁺/Cl^- ionic/covalent curve crossing ============================================================================

Status:: Root cause identified, not yet fixed
Found in:: Job 3867 on the dev dashboard – Na⁺…Cl^- 40-point log-spaced R-scan (1.30-20.00 A), REVDSD-PBEP86-D4/2021 def2-TZVPPD, --charges "1 -1", through the M3 ORCA N-fragment substep this campaign built. Structures: ~/Sites/training-data/2026-08-03_na-cl-scan-run1/ (Na-Cl.sdf).
Affects:: The cluster sub-calculation only (the one full-complex, both atoms real, net-neutral SCF – one of the 2N + 1 jobs per point). Not the fragment-alone or fragment-in-cluster-basis jobs.
Reported by:: user, from the labeled Dimers.extxyz – noticed nonzero x/y forces that should vanish by the molecule’s linear symmetry, plus noisy Chargemol/DDEC6 charges with no clean trend vs. R.

Summary#

This is not noise and not a convergence-tolerance problem. A single- determinant (RHF/KS) description of neutral NaCl has two nearly-degenerate electronic solutions once R is large enough – an ionic one (Na⁺Cl^-, the desired diabatic surface for this training set) and a covalent one (Na^.Cl^., two neutral atoms). They cross where the Coulomb stabilization of the ionic state equals the cost of forming it:

\[R_c \approx \frac{e^2}{\mathrm{IP}(\mathrm{Na}) - \mathrm{EA}(\mathrm{Cl})} = \frac{14.399\ \mathrm{eV \cdot A}}{5.139 - 3.617\ \mathrm{eV}} \approx 9.5\ \text{A}\]

which lands right in the middle of the observed problem window. Near a crossing like this, the SCF’s default (neutral-atom) initial guess can converge – self-consistently, no convergence warnings – to either branch depending on details of the specific geometry; nothing chains guesses between R-scan points (checked: no %moinp/MORead in any orca.inp, each point is a fresh default guess), so which branch wins is essentially decided independently, geometry by geometry. Tighter SCF thresholds do not fix this: the SCF is converging, just to the wrong one of two locally-stable solutions.

Evidence#

Pulled the per-subcalc ORCA logs for job 3867 (root/5/<config>/2/1/ {cluster,1-alone,2-alone,1-in-cluster,2-in-cluster}/orca.out) for the flagged points and several clean neighbors for comparison.

Isolated to the cluster job. The “cluster” SCF energy is non-monotonic in R = 6.5-12.2 A (re-deepens instead of relaxing toward the E(Na⁺) + E(Cl^-) asymptote), then snaps back to sane behavior at R = 13.13 A. Over the identical R range, the 1-alone, 2-alone, 1-in-cluster, and 2-in-cluster energies are all smooth to ~1e-6 Hartree – unsurprising, since each of those is a single, unambiguous charge state (+1, -1, or a real fragment plus inert ghost functions with no competing electron count), so there is no second solution for them to fall into.
A virtual orbital is dropping toward becoming an electron-transfer acceptor state. In the cluster job, the LUMO(+1) pair’s orbital energy drops steadily from +0.0147 Ha at R = 6.52 A through negative values by R = 13.13 A – the textbook signature of approaching this kind of crossing.
Genuine SCF symmetry breaking. The molecule is linear (D:math:_{infty h}), so degenerate pi-orbital pairs must be exactly degenerate. At R = 12.25 A (config 1,33) the normally-exact pair splits slightly (-5.4467 vs. -5.4436 eV) – the same phenomenon that shows up downstream as nonzero x/y forces.
Total-force check pinpoints the exact broken points. Sum of forces over the two atoms must be exactly zero (translational invariance, no external field); this is a stronger, purely computational diagnostic than “x/y should vanish by symmetry” and needs no chemistry knowledge to apply. Background level (healthy points) is ~1e-4 to 1e-3; three points are 50-500x above that:

config

R (A)

|sum F|

in the crossing window?

1,17

3.99

0.049

no – likely an unrelated, isolated glitch

1,25

6.99

0.066

yes (edge)

1,27

8.04

0.054

yes

The energy is smooth through all three points – only the gradient (and, for 1,25/1,27, the Chargemol charge, and the interaction energy itself in Results.csv) is affected. The broader charge/energy drift spans roughly configs 1,26-1,33 (R = 7.5-12.25 A) even though only two of those points trip the force-sum threshold – i.e. the crossing can perturb a point’s density/charge partition measurably before it perturbs the gradient enough to trip a simple threshold, so a threshold-based gate alone will under-count affected points. 1,17 at R = 3.99 A is well inside the normal ionic-bonding region, nowhere near \(R_c\) – treat it as a separate, one-off SCF quirk rather than the same phenomenon until shown otherwise.

config	R (A)	\|sum F\|	in the crossing window?
`1,17`	3.99	0.049	no – likely an unrelated, isolated glitch
`1,25`	6.99	0.066	yes (edge)
`1,27`	8.04	0.054	yes

Why this matters beyond Na/Cl#

\(R_c\) depends only on IP(cation) - EA(anion), so every ionic pair this campaign scans will have some crossing distance, just at a different R. Worth precomputing \(R_c\) (or at least flagging “this pair’s IP-EA gap is small enough that \(R_c\) falls inside the planned scan range”) before each new ion-pair production run, rather than rediscovering it after the fact per-pair. Ion-water and ion-water-water systems (the angular/trimer test sets built alongside this scan, see ~/Sites/mlff-training/2026-08-07_ion-water-angular-and-trimers/) don’t have this exact failure mode since the “cluster” fragment count/composition differs, but the same class of problem – a fragment sub-calculation with more than one accessible electronic state near-degenerate at some geometry – is worth keeping in mind generally for charged multi-fragment BSSE.

Fix directions (not yet implemented)#

Because this is “converges to the wrong stable solution,” not “fails to converge,” a tighter TightSCF -> VeryTightSCF will not help. In order of expected leverage:

Seed the cluster job’s guess from the fragment orbitals already being computed. Every BSSE point already runs 1-alone (Na⁺) and 2-alone (Cl^-) as two of the 2N + 1 required jobs – their orbitals are unambiguous (single charge state each) and free. Combine them into the cluster job’s initial guess (ORCA supports building a guess from fragment MOs) instead of letting it fall back to the default neutral-atom guess. This directly biases every point toward the ionic diabatic surface, which is already the deliberate target everywhere else in this pipeline (the whole point of --charges fixing per-fragment charge), including past \(R_c\) where the true adiabatic ground state would actually be covalent – that’s correct for training- set purposes, not a bug to route around.
Add SCF stability analysis (!StabPerform in current ORCA syntax) to the cluster job as a backstop, so a solution that is actually a saddle point gets kicked back toward a genuine minimum. Cheaper to add than (1) but doesn’t by itself guarantee landing on the ionic minimum rather than the covalent one past \(R_c\) – (1) and (2) are complementary, not alternatives.
A generic post-hoc QC gate, independent of the above: check \(\sum \vec F \approx 0\) across the full cluster automatically for every point of every charged/multi-fragment production run (cheap, no chemistry-specific knowledge needed) and flag/retry anything outside tolerance. This is what surfaced the problem here and generalizes to any future ionic system, but per the point above about 1,26/1,28-1,33 not all tripping the threshold, it should be a first-pass filter, not the only check – pair it with a smoothness check on the CP-corrected energy and/or the atomic charges vs. neighboring R points.

Raw data#

Dashboard: dev (http://localhost:55066), job 3867, project water. Key files: Dimers.extxyz (labeled structures/forces/charges), Results.csv (per-point interaction energies), and per-subcalc ORCA I/O under root/5/<config>/2/1/{cluster,1-alone,2-alone,1-in-cluster, 2-in-cluster}/orca.{inp,out} for each <config> (e.g. 1,25).