Add SCC convergence diagnostics and benchmark

Author: QG-phy

.gitignore

@@ -46,6 +46,7 @@ __pycache__/
 *.py[cod]
 *$py.class
 .vscode
+.codegraph/
 # C extensions
 *.so
 dptb/tests/data/esk_orb/sktb.json

docs/maintenance/scc_numerics_test_plan.md

@@ -0,0 +1,210 @@
+# Carbon BCC Self-Consistent Total-Energy K-Point Convergence Test Plan
+
+## Status
+
+In progress. The current scope is the convergence of the self-consistent total
+energy of one standalone Carbon BCC calculation with respect to k-point
+density. Inner SCC charge convergence is a required diagnostic but is not the
+target convergence criterion. This is a validation plan, not an implementation
+proposal.
+
+## Motivation
+
+DeePTB's SCC-DFTB functionality should give a stable, reproducible total energy
+for a fixed Carbon BCC input as the k-point density is increased. An individual
+mesh reaching its SCC charge tolerance only proves that mesh was solved
+self-consistently; it does not prove Brillouin-zone integration or total-energy
+convergence. This work is not a direct test against the existing multi-
+structure DFTB+ pytest benchmark.
+
+The downstream UniSK workflow exposed a practical concern. Its legacy Carbon
+Bayesian-optimization script uses per-structure k-point spacings such as
+`[0.05, 0.05, 0.05]` for BCC and FCC. For the Carbon BCC cell used by the
+DeePTB DFTB+ benchmark, that spacing generates a `75 x 75 x 75` mesh
+(`421875` raw points), whereas the existing DeePTB benchmark uses
+`20 x 20 x 20` (`8000` raw points). The former was copied from an old
+production-oriented workflow with the comment that it was chosen to "ensure
+convergence"; no k-point convergence study was preserved with that script.
+
+This difference must not be interpreted as proof that DeePTB SCC is slow or
+incorrect. It is first a numerical-configuration question. A mesh that is
+overly fine can make a correct SCC calculation impractical inside a Bayesian
+optimization loop.
+
+## Existing Assets
+
+The repository contains a broader correctness regression:
+
+- `dptb/tests/test_dftb_scc.py::test_dftbscc_matches_dftbp_benchmarks`
+  compares DeePTB SCC electronic energies against DFTB+ reference data for
+  Carbon BCC, graphene, and dimer EOS scans.
+- The scan contains seven lattice/bond-length scales for each structure.
+- DFTB+ reference tables and the corresponding structures/SK files live in
+  `dptb/tests/data/dftb/structs_eos` and `dptb/tests/data/dftb`.
+- The current acceptance threshold is a maximum electronic-energy difference
+  below `1e-2 eV` at every scale.
+- Its direct command is:
+
+  ```bash
+  python -m pytest -q dptb/tests/test_dftb_scc.py -k matches_dftbp_benchmarks
+  ```
+
+That regression is useful background and supplies the canonical Carbon
+structure and SK data, but it is not the execution target for the present
+study. The current standalone entry point is:
+
+- `examples/dftb_scc/carbon_bcc_scc_convergence/run.py`;
+- `examples/dftb_scc/carbon_bcc_scc_convergence/POSCAR`, containing only the
+  scale-1.000 primitive Carbon BCC cell; and
+- the canonical `C-C.skf` under `dptb/tests/data/dftb`, shared without invoking
+  any pytest helper or multi-structure EOS scan.
+
+The standalone example records self-consistent total, electronic, and repulsive
+energies; energy changes between successive meshes and against all denser
+meshes in the scan; k-point counts; SCC iteration count and charge residual
+history; electron count; Mulliken charge; Fermi level; and wall time. Existing
+unit tests for SCC state, mixers, API equivalence, and k-point symmetry
+reduction remain useful but do not define this study's convergence criterion.
+
+## Questions To Answer
+
+The following effects are coupled in an SCC calculation but need separate
+tests and diagnostics:
+
+1. **K-point integration:** What mesh keeps the self-consistent total energy per
+   atom within the agreed tolerance of every denser mesh tested for the fixed
+   Carbon BCC cell?
+2. **Occupations and smearing:** Are Fermi-Dirac and Gaussian paths correctly
+   normalized, electron-number conserving, and stable in the relevant
+   zero-/finite-temperature limits?
+3. **Fermi-level solve:** Does the root solve use the same k-point weights and
+   occupation convention as the band-energy and charge calculations?
+4. **Symmetry reduction:** Do full, time-reversal-reduced, and
+   rotational-symmetry-reduced meshes give equivalent observables with correct
+   accumulated weights?
+5. **Inner SCC loop:** Does every k mesh reach the charge-residual tolerance,
+   and can SCC numerical noise contaminate the outer total-energy trend?
+6. **Performance:** Which component dominates runtime after the numerical
+   configuration is fixed: Hamiltonian construction, eigensolve, occupation,
+   Mulliken population, Gamma/SCC shift, or mixing?
+
+## Important Configuration Distinction
+
+The current UniSK Carbon workflow and the DeePTB DFTB+ benchmark are not a
+direct apples-to-apples timing or numerical comparison. Among other settings,
+the legacy UniSK path uses Gaussian smearing (`"G"`) and a `0.05` spacing,
+while the DFTB+ benchmark uses Fermi-Dirac smearing (`"FD"`) and fixed meshes
+such as Carbon BCC `20 x 20 x 20`. Future comparisons must explicitly record
+all physical and numerical settings, including SK files, repulsive parameters,
+cutoffs, temperature, mesh, symmetry flags, mixer, tolerance, and maximum
+iteration count.
+
+## Phased Work
+
+### Phase 0: Lock the standalone Carbon BCC baseline
+
+- Use only the scale-1.000 Carbon BCC structure in the standalone example.
+- Do not call `test_dftbscc_matches_dftbp_benchmarks` or scan graphene, dimer,
+  or the seven-point Carbon EOS in this phase.
+- Record wall time, SCC iteration count, full residual history, final residual,
+  electron count, Fermi level, total/electronic/repulsive energies, SCC shift
+  energy, and Mulliken charge for every numerical case.
+- Keep the command runnable independently of notebooks and without generating
+  or duplicating SK tables.
+
+### Phase 1: K-point convergence matrix
+
+For the fixed Carbon BCC cell, scan deliberately chosen fixed meshes and, where
+useful, equivalent spacing-derived meshes. At minimum record:
+
+- total-energy change per atom from the preceding mesh;
+- maximum absolute total-energy change per atom relative to every denser mesh
+  in the same run;
+- an explicit k-point convergence flag evaluated against a documented energy
+  threshold; the last and densest point cannot validate itself;
+- charge and Fermi-level change relative to the finest accepted mesh;
+- SCC iteration count and convergence failures;
+- raw and symmetry-reduced k-point counts;
+- wall time split by SCC iteration and total calculation.
+
+Use the results to define separate **smoke**, **development**, and
+**production** mesh presets. Do not retain `0.05` as a default solely because
+it appeared in a legacy script.
+
+The initial `20^3, 24^3, 28^3, 32^3, 36^3, 40^3` even-grid scan shows that
+`20^3` differs from a denser result by up to `2.252 meV/atom`. It therefore does
+not meet a `1 meV/atom` criterion. `24^3` and denser tested even grids remain
+within `1 meV/atom`, but neighbouring odd grids and spacing-derived meshes must
+still be checked before fixing a production preset.
+
+### Phase 2: Occupation and Fermi-level regression tests
+
+- Add small, deterministic tests for occupation normalization and exact
+  electron-number conservation with non-uniform k-point weights.
+- Test Fermi-Dirac and Gaussian implementations independently.
+- Verify expected low-temperature and high-temperature limiting behaviour.
+- First record Fermi level, band energy, SCC charge, and electron-count error
+  for the standalone Carbon BCC input. Broader material-class coverage is a
+  separate follow-up.
+
+### Phase 3: Symmetry-reduction equivalence
+
+- For the Carbon BCC cell, compare full meshes with time-reversal and
+  rotational reductions.
+- Assert equivalence of electron count, Fermi level, charge, total/electronic
+  energy, SCC shift energy, and converged SCC shift within explicit tolerances.
+- Assert that the sum of reduced k-point weights equals the full mesh weight.
+- Report the reduction factor as a performance metric.
+
+Symmetry reduction is a performance optimization, not a change to the target
+physical result.
+
+### Phase 4: SCC robustness and performance work
+
+Only after Phases 0-3 identify the limiting configuration should solver work
+begin. Candidate changes include static-data caching across SCC iterations,
+batched k-point eigensolves, warm starts, and mixer/preconditioner tuning.
+Every change must run the numerical regression matrix and a separate timing
+benchmark.
+
+The primitive Carbon BCC cell has one atom, so charge conservation and symmetry
+largely pin its Mulliken charge. It can reveal occupation-noise effects in a
+very tight SCC tolerance, but it is not a difficult charge-mixing problem. A
+true mixer stress test will require a separate cell with non-equivalent charge
+degrees of freedom or an explicitly perturbed initial charge; that must not be
+silently substituted into this Carbon BCC baseline.
+
+## Tetrahedron Integration
+
+Tetrahedron integration may be useful for zero-temperature Brillouin-zone
+integration, DOS, and difficult metallic cases, but it is not a drop-in
+replacement for smearing. It requires a separate design covering band
+interpolation, degeneracies, weights, self-consistency, and eventually force
+consistency. It should be evaluated only after the existing Fermi-Dirac and
+Gaussian paths are fully benchmarked and tested.
+
+## Acceptance Criteria For Future Changes
+
+A proposed SCC numerical or performance change is acceptable only when:
+
+1. the standalone Carbon BCC total-energy matrix remains within its documented
+   meV/atom tolerance against all denser meshes tested, while inner SCC charge
+   residuals also meet their separate tolerance;
+2. symmetry-reduced and full-mesh results agree within documented tolerances;
+3. electron-number conservation is verified for every occupation path;
+4. convergence failures and iteration counts do not regress without a recorded
+   numerical reason; and
+5. any claimed speedup reports the mesh, symmetry flags, precision, hardware,
+   warm-up policy, and the same numerical acceptance criteria.
+
+## Non-Goals
+
+- Covering graphene, the Carbon dimer, or the seven-point Carbon EOS in the
+  current standalone SCC convergence study.
+- Using `test_dftbscc_matches_dftbp_benchmarks` as the current execution entry
+  point.
+- Changing default UniSK production parameters in this task.
+- Declaring a legacy `kmeshspacing` value physically necessary without a
+  convergence study.
+- Replacing the SCC algorithm with a tetrahedron method before the current
+  occupation and k-point paths have a complete regression suite.

dptb/nn/dftb/dftb_scc.py

@@ -85,7 +85,7 @@ def reset(self) -> None:
         structures, as it avoids the overhead of reinitializing the expensive
         model, SK parameters, and transformers for each structure.
 
-        The following attributes are reset to None:
+        The following per-calculation state is cleared or reset:
         - atomic_numbers: Atomic numbers of the current structure
         - elec_totE, elec_H0_bandE, elec_bandE: Electronic energies
         - E_fermi: Fermi energy
@@ -96,6 +96,9 @@ def reset(self) -> None:
         - inv_r: Inverse distance matrix
         - Gamma: Coulomb interaction matrix
         - total_energy, total_rep_energy: Total energies
+        - n_iterations: Number of Mulliken/SCC iterations performed
+        - charge_residual_history: Maximum charge residual from each SCC update
+        - final_charge_residual: Last maximum charge residual, if one was evaluated
 
         Also resets the Mulliken calculator state via mulliken.reset().
         '''
@@ -115,6 +118,9 @@ def reset(self) -> None:
         self.total_energy = None
         self.total_rep_energy = None
         self.is_converged = False  # Track SCC convergence status
+        self.n_iterations = 0
+        self.charge_residual_history = []
+        self.final_charge_residual = None
         # Also reset the Mulliken calculator state
         self.mulliken.reset()
 
@@ -218,6 +224,12 @@ def run_iters(self,
             - self.elec_totE: Total electronic energy (eV)
             - self.E_fermi: Fermi energy (eV)
             - self.mulliken.mul_charge: Final Mulliken charges
+            - self.is_converged: Whether the charge residual reached ``tol``
+            - self.n_iterations: Number of Mulliken/SCC iterations performed
+            - self.charge_residual_history: Maximum charge residual for each
+              update after the initial charge evaluation
+            - self.final_charge_residual: Last value in the residual history,
+              or ``None`` when no residual was evaluated
 
         '''
 
@@ -328,10 +340,13 @@ def run_iters(self,
 
                 # Convergence criterion: ||q_diff|| = ||q_out - q_inp|| (DFTBplus convention)
                 diff = np.abs(q_diff)
-                log.debug(f' ITERATION: {int(iteration)}  Max charge difference: {diff.max().item():.15f} e')
+                max_charge_residual = float(diff.max())
+                self.charge_residual_history.append(max_charge_residual)
+                self.final_charge_residual = max_charge_residual
+                log.debug(f' ITERATION: {int(iteration)}  Max charge difference: {max_charge_residual:.15f} e')
                 iteration += 1
 
-                if diff.max() < tol: # Converged and not exiting due to max_iter
+                if max_charge_residual < tol: # Converged and not exiting due to max_iter
                     assert self.scc_shift_energy is not None, "scc_shift_energy should not be None when converged."
                     assert self.elec_H0_bandE is not None, "elec_H0_bandE should not be None when converged."
                     log.debug(f'Convergence reached after {int(iteration)} iters.')
@@ -374,6 +389,7 @@ def run_iters(self,
             data[AtomicDataDict.EIGENVECTOR_KEY] = torch.cat(eigvecs, dim=0)
 
         self.data = data
+        self.n_iterations = iteration
         self.E_fermi = self.mulliken.estimated_E_fermi
         self.elec_bandE = self.mulliken.elec_bandE
         self.elec_totE = self.elec_H0_bandE + self.scc_shift_energy

dptb/tests/test_dftb_scc.py

@@ -478,6 +478,9 @@ def test_dftbscc_reset_clears_state(rootdir=rootdir):
     assert dftbscc.expGamma_onsite is not None, "expGamma_onsite should be set after calculation"
     assert dftbscc.inv_r is not None, "inv_r should be set after calculation"
     assert dftbscc.Gamma is not None, "Gamma should be set after calculation"
+    assert dftbscc.n_iterations > 0, "n_iterations should be populated after calculation"
+    assert len(dftbscc.charge_residual_history) == dftbscc.n_iterations - 1
+    assert dftbscc.final_charge_residual == dftbscc.charge_residual_history[-1]
 
     # Call reset explicitly
     dftbscc.reset()
@@ -498,6 +501,9 @@ def test_dftbscc_reset_clears_state(rootdir=rootdir):
     assert dftbscc.Gamma is None, "Gamma should be None after reset"
     assert dftbscc.total_energy is None, "total_energy should be None after reset"
     assert dftbscc.total_rep_energy is None, "total_rep_energy should be None after reset"
+    assert dftbscc.n_iterations == 0, "n_iterations should be zero after reset"
+    assert dftbscc.charge_residual_history == [], "residual history should be empty after reset"
+    assert dftbscc.final_charge_residual is None, "final residual should be None after reset"
 
     # Verify persistent state (model, skp, etc.) is NOT cleared
     assert dftbscc.model is not None, "model should persist after reset"

examples/dftb_scc/README.md

@@ -6,6 +6,10 @@ Imports and obvious absolute/local SK paths were adjusted to use DeePTB modules
 
 Recommended notebooks:
 
+- `carbon_bcc_scc_convergence/run.py`: standalone Carbon BCC self-consistent
+  total-energy convergence against k meshes/spacings. Inner SCC, mixer, and
+  smearing diagnostics are recorded separately. See the example README for
+  commands and interpretation.
 - `benchmark_eos_dftbp.ipynb`: SCC EOS comparison against DFTB+ reference data in `dptb/tests/data/dftb`.
 - `hBN_scc/test_scc_hBN.ipynb`: hBN SCC workflow and band calculation.
 - `hBN_scc/test_scc_hBN_nnsk.ipynb`: hBN nonorthogonal NNSK-SCC workflow.

examples/dftb_scc/carbon_bcc_scc_convergence/POSCAR

@@ -0,0 +1,9 @@
+C BCC primitive cell, scale 1.000
+1.0000000000000000
+  2.0557030981993387   0.0000005029268204   0.0000008733164838
+ -0.6852331938919997   1.9381358756482154   0.0000008681559103
+ -0.6852331920572658  -0.9690657046357299   1.6784761976400846
+C
+1
+Cartesian
+ -0.0000000006781663   0.0000000005972072  -0.0000000000918888