Correct SCF-step convergence residuals and split delta_density_rms (#452, #453, #454)

Fix `delta_energies_total` to derive from the per-SCF `scf_steps.energies_total`
series (not the labeled `Outputs.total_energies`), and stop synthesizing
`delta_force_abs` from the final `total_forces` (now parser-only).

Replace `delta_density_rms` with definition-specific quantities -- `delta_charge_abs`,
`delta_charge_density_rms`, `delta_charge_relative`, `delta_density_matrix_rms`,
`delta_density_matrix_max` -- plus `energy_error_estimate`, and replace
`ChargeConvergenceTarget` with a `DensityConvergenceTarget` whose `type` selector
derives the path and unit from the referenced `scf_steps` field.

Breaking: removes `delta_density_rms`; renames `ChargeConvergenceTarget`.

Author: ndaelman-hu

docs/schema/index.md

@@ -170,7 +170,7 @@ Thermodynamics workflow for free-energy and thermodynamic property calculations
 
 Convergence target classes and workflow-level convergence result structures
 
-**Key sections:** WorkflowConvergenceTarget, EnergyConvergenceTarget, ForceConvergenceTarget, PotentialConvergenceTarget, ChargeConvergenceTarget, WavefunctionConvergenceTarget, WorkflowConvergenceResults, SimulationWorkflowModel, SimulationWorkflowResults, GeometryOptimizationModel, GeometryOptimizationResults
+**Key sections:** WorkflowConvergenceTarget, EnergyConvergenceTarget, ForceConvergenceTarget, PotentialConvergenceTarget, DensityConvergenceTarget, WavefunctionConvergenceTarget, WorkflowConvergenceResults, SimulationWorkflowModel, SimulationWorkflowResults, GeometryOptimizationModel, GeometryOptimizationResults
 
 ## [Workflow Trajectory Properties](workflow_trajectory.md)
 

docs/schema/outputs.diagram.md

@@ -15,44 +15,44 @@ This diagram shows the relationships between schema classes:
 classDiagram
     class AbsorptionSpectrum {
     }
-    class ChemicalPotential {
-    }
-    class CrystalFieldSplitting {
+    class ElectronicBandGap {
     }
-    class ElectronicBandStructure {
+    class ElectronicEigenvalues {
     }
-    class ElectronicDensityOfStates {
+    class ElectronicGreensFunction {
     }
     class ElectronicSelfEnergy {
     }
-    class FermiSurface {
+    class HoppingMatrix {
+    }
+    class KineticEnergy {
     }
     class Occupancy {
     }
     class Outputs {
     }
     class PhysicalProperty {
     }
-    class PotentialEnergy {
+    class SCFSteps {
     }
-    class QuasiparticleWeight {
+    class Temperature {
     }
-    class TotalForce {
+    class TotalEnergy {
     }
-    class XASSpectrum {
+    class TotalForce {
     }
     Outputs *-- AbsorptionSpectrum : absorption_spectra
-    Outputs *-- ChemicalPotential
-    Outputs *-- CrystalFieldSplitting
-    Outputs *-- ElectronicBandStructure
-    Outputs *-- ElectronicDensityOfStates : electronic_dos
+    Outputs *-- ElectronicBandGap
+    Outputs *-- ElectronicEigenvalues
+    Outputs *-- ElectronicGreensFunction
     Outputs *-- ElectronicSelfEnergy : electronic_self_energies
-    Outputs *-- FermiSurface
+    Outputs *-- HoppingMatrix : hopping_matrices
+    Outputs *-- KineticEnergy : kinetic_energies
     Outputs *-- Occupancy : occupancies
-    Outputs *-- PotentialEnergy : potential_energies
-    Outputs *-- QuasiparticleWeight
+    Outputs *-- SCFSteps
+    Outputs *-- Temperature
+    Outputs *-- TotalEnergy : total_energies
     Outputs *-- TotalForce
-    Outputs *-- XASSpectrum : xas_spectra
 ```
 
 </div>
@@ -67,43 +67,43 @@ _Diagram 2 of 2 (split due to large number of children)_
 
 ```mermaid
 classDiagram
-    class ElectronicBandGap {
+    class ChemicalPotential {
     }
-    class ElectronicEigenvalues {
+    class CrystalFieldSplitting {
     }
-    class ElectronicGreensFunction {
+    class ElectronicBandStructure {
     }
-    class HoppingMatrix {
+    class ElectronicDensityOfStates {
     }
-    class HybridizationFunction {
+    class FermiSurface {
     }
-    class KineticEnergy {
+    class HybridizationFunction {
     }
     class Outputs {
     }
     class Permittivity {
     }
     class PhysicalProperty {
     }
-    class RadiusOfGyration {
+    class PotentialEnergy {
     }
-    class SCFSteps {
+    class QuasiparticleWeight {
     }
-    class Temperature {
+    class RadiusOfGyration {
     }
-    class TotalEnergy {
+    class XASSpectrum {
     }
-    Outputs *-- ElectronicBandGap
-    Outputs *-- ElectronicEigenvalues
-    Outputs *-- ElectronicGreensFunction
-    Outputs *-- HoppingMatrix : hopping_matrices
+    Outputs *-- ChemicalPotential
+    Outputs *-- CrystalFieldSplitting
+    Outputs *-- ElectronicBandStructure
+    Outputs *-- ElectronicDensityOfStates : electronic_dos
+    Outputs *-- FermiSurface
     Outputs *-- HybridizationFunction
-    Outputs *-- KineticEnergy : kinetic_energies
     Outputs *-- Permittivity : permittivities
+    Outputs *-- PotentialEnergy : potential_energies
+    Outputs *-- QuasiparticleWeight
     Outputs *-- RadiusOfGyration : radii_of_gyration
-    Outputs *-- SCFSteps
-    Outputs *-- Temperature
-    Outputs *-- TotalEnergy : total_energies
+    Outputs *-- XASSpectrum : xas_spectra
 ```
 
 <p class="uml-legend__title">Legend</p>

docs/schema/outputs.md

@@ -91,8 +91,13 @@ classDiagram
 |---|---|---|
 | `energies_total` | m_float64(float) (shape: ['*']) | Total energy at each SCF step. |
 | `delta_energies_total` | m_float64(float) (shape: ['*']) | Absolute change of total energy at each SCF step. |
+| `energy_error_estimate` | m_float64(float) (shape: ['*']) | <details><summary>Estimate of the remaining error in the total energy at each SCF step,</summary>Estimate of the remaining error in the total energy at each SCF step,<br>derived from the density residual rather than from the change of the<br>total energy itself. For example, Quantum ESPRESSO's "estimated scf<br>accuracy" is the Hartree self-energy of the density residual. Distinct<br>from `delta_energies_total`, which is the change of the total energy<br>between consecutive steps.</details> |
 | `delta_potential_rms` | m_float64(float) (shape: ['*']) | Root mean square of change of potential energy at each SCF step. |
-| `delta_density_rms` | m_float64(float) (shape: ['*']) | Root mean square of change of potential energy at each SCF step. |
+| `delta_charge_abs` | m_float64(float) (shape: ['*']) | <details><summary>Volume-integrated absolute change of the electron density between</summary>Volume-integrated absolute change of the electron density between<br>consecutive SCF steps, `integral \|rho_n(r) - rho_(n-1)(r)\| d^3r`,<br>expressed as a charge (equivalently a number of electrons). Reported by<br>all-electron codes such as WIEN2k (`:DIS`). The exact norm and any<br>normalization are a code-reported convention that the schema does not<br>enforce.</details> |
+| `delta_charge_density_rms` | m_float64(float) (shape: ['*']) | <details><summary>Root mean square, over real-space grid points, of the change of the</summary>Root mean square, over real-space grid points, of the change of the<br>electron density between consecutive SCF steps. Unlike `delta_charge_abs`<br>the volume is retained, so this is a charge density. Reported by<br>plane-wave codes such as VASP (`rms(c)`).</details> |
+| `delta_charge_relative` | m_float64(float) (shape: ['*']) | Integrated absolute density change normalized by the electron count, `integral \|rho_n - rho_(n-1)\| d^3r / N`, hence dimensionless. Reported by exciting ("charge distance") and GPAW (per valence electron). |
+| `delta_density_matrix_rms` | m_float64(float) (shape: ['*']) | <details><summary>Root mean square of the change of the density-matrix elements `P_munu`</summary>Root mean square of the change of the density-matrix elements `P_munu`<br>(in the non-orthonormal atomic-orbital basis) between consecutive SCF<br>steps. The elements are dimensionless, so is this residual. Reported by<br>Gaussian-basis codes such as CRYSTAL (`tst`) and ORCA (`RMS-DP`).</details> |
+| `delta_density_matrix_max` | m_float64(float) (shape: ['*']) | <details><summary>Maximum absolute change of the density-matrix elements `P_munu` between</summary>Maximum absolute change of the density-matrix elements `P_munu` between<br>consecutive SCF steps; the max-norm counterpart of<br>`delta_density_matrix_rms`. Reported by Gaussian-basis codes such as<br>CP2K, CRYSTAL (`PX`), and ORCA (`Max-DP`).</details> |
 | `delta_wavefunction_rms` | m_float64(float) (shape: ['*']) | Root mean square of change of wavefunction coefficients at each SCF step. Dimensionless quantity representing convergence of orbital coefficients. |
 | `delta_force_abs` | m_float64(float) (shape: ['*']) | Absolute change of forces at each SCF step. |
 | `durations` | m_float64(float) (shape: ['*']) | Time spent at each SCF step. |

docs/schema/workflow_convergence.diagram.md

@@ -13,7 +13,7 @@ This diagram shows the relationships between schema classes:
 
 ```mermaid
 classDiagram
-    class ChargeConvergenceTarget {
+    class DensityConvergenceTarget {
     }
     class EnergyConvergenceTarget {
     }
@@ -35,7 +35,7 @@ classDiagram
     }
     class WorkflowConvergenceTarget {
     }
-    WorkflowConvergenceTarget <|-- ChargeConvergenceTarget
+    WorkflowConvergenceTarget <|-- DensityConvergenceTarget
     WorkflowConvergenceTarget <|-- EnergyConvergenceTarget
     WorkflowConvergenceTarget <|-- ForceConvergenceTarget
     SimulationWorkflowResults <|-- GeometryOptimizationResults

docs/schema/workflow_convergence.md

@@ -10,7 +10,7 @@
 
 ```mermaid
 classDiagram
-    class ChargeConvergenceTarget
+    class DensityConvergenceTarget
     class EnergyConvergenceTarget
     class ForceConvergenceTarget
     class GeometryOptimizationModel
@@ -21,7 +21,7 @@ classDiagram
     class WavefunctionConvergenceTarget
     class WorkflowConvergenceResults
     class WorkflowConvergenceTarget
-    WorkflowConvergenceTarget <|-- ChargeConvergenceTarget
+    WorkflowConvergenceTarget <|-- DensityConvergenceTarget
     WorkflowConvergenceTarget <|-- EnergyConvergenceTarget
     WorkflowConvergenceTarget <|-- ForceConvergenceTarget
     SimulationWorkflowResults <|-- GeometryOptimizationResults
@@ -84,15 +84,16 @@ classDiagram
 |---|---|---|
 | `threshold` | m_float_bounded(float) | <details><summary>Convergence threshold.</summary>Convergence threshold. Must be non-negative.<br>When threshold_type is 'relative', must be dimensionless.<br>When threshold_type is 'absolute', 'maximum', or 'rms', must have physical units.<br>Child classes override this to add convergence path annotations.</details> |
 
-### `ChargeConvergenceTarget`
+### `DensityConvergenceTarget`
 
 | Section | Description | MetaInfo |
 |---|---|---|
-| `ChargeConvergenceTarget` | Convergence target for electron density/charge differences. | [Open in MetaInfo browser](https://nomad-lab.eu/prod/v1/develop/gui/analyze/metainfo/nomad_simulations/section_definitions@nomad_simulations.schema_packages.workflow.general.ChargeConvergenceTarget){:target="_blank"} |
+| `DensityConvergenceTarget` | Convergence target for the SCF electron-density residual. | [Open in MetaInfo browser](https://nomad-lab.eu/prod/v1/develop/gui/analyze/metainfo/nomad_simulations/section_definitions@nomad_simulations.schema_packages.workflow.general.DensityConvergenceTarget){:target="_blank"} |
 
 | Quantity | Type | Description |
 |---|---|---|
 | `threshold` | m_float_bounded(float) | <details><summary>Convergence threshold.</summary>Convergence threshold. Must be non-negative.<br>When threshold_type is 'relative', must be dimensionless.<br>When threshold_type is 'absolute', 'maximum', or 'rms', must have physical units.<br>Child classes override this to add convergence path annotations.</details> |
+| `type` | Enum | Which density-convergence residual this target tracks. Selects the `scf_steps.delta_<type>` quantity read for the check (and thereby its expected unit). |
 
 ### `WavefunctionConvergenceTarget`
 

scripts/verticals.py

@@ -70,7 +70,7 @@
             'EnergyConvergenceTarget',
             'ForceConvergenceTarget',
             'PotentialConvergenceTarget',
-            'ChargeConvergenceTarget',
+            'DensityConvergenceTarget',
             'WavefunctionConvergenceTarget',
             'WorkflowConvergenceResults',
             'SimulationWorkflowModel',

src/nomad_simulations/schema_packages/outputs.py

@@ -64,6 +64,20 @@ class SCFSteps(ArchiveSection):
         """,
     )
 
+    energy_error_estimate = Quantity(
+        shape=['*'],
+        type=float,
+        unit='joule',
+        description="""
+        Estimate of the remaining error in the total energy at each SCF step,
+        derived from the density residual rather than from the change of the
+        total energy itself. For example, Quantum ESPRESSO's "estimated scf
+        accuracy" is the Hartree self-energy of the density residual. Distinct
+        from `delta_energies_total`, which is the change of the total energy
+        between consecutive steps.
+        """,
+    )
+
     delta_potential_rms = Quantity(
         shape=['*'],
         type=float,
@@ -73,12 +87,64 @@ class SCFSteps(ArchiveSection):
         """,
     )
 
-    delta_density_rms = Quantity(
+    delta_charge_abs = Quantity(
         shape=['*'],
         type=float,
         unit='coulomb',
         description="""
-        Root mean square of change of potential energy at each SCF step.
+        Volume-integrated absolute change of the electron density between
+        consecutive SCF steps, `integral |rho_n(r) - rho_(n-1)(r)| d^3r`,
+        expressed as a charge (equivalently a number of electrons). Reported by
+        all-electron codes such as WIEN2k (`:DIS`). The exact norm and any
+        normalization are a code-reported convention that the schema does not
+        enforce.
+        """,
+    )
+
+    delta_charge_density_rms = Quantity(
+        shape=['*'],
+        type=float,
+        unit='coulomb / meter ** 3',
+        description="""
+        Root mean square, over real-space grid points, of the change of the
+        electron density between consecutive SCF steps. Unlike `delta_charge_abs`
+        the volume is retained, so this is a charge density. Reported by
+        plane-wave codes such as VASP (`rms(c)`).
+        """,
+    )
+
+    delta_charge_relative = Quantity(
+        shape=['*'],
+        type=float,
+        unit='dimensionless',
+        description="""
+        Integrated absolute density change normalized by the electron count,
+        `integral |rho_n - rho_(n-1)| d^3r / N`, hence dimensionless. Reported by
+        exciting ("charge distance") and GPAW (per valence electron).
+        """,
+    )
+
+    delta_density_matrix_rms = Quantity(
+        shape=['*'],
+        type=float,
+        unit='dimensionless',
+        description="""
+        Root mean square of the change of the density-matrix elements `P_munu`
+        (in the non-orthonormal atomic-orbital basis) between consecutive SCF
+        steps. The elements are dimensionless, so is this residual. Reported by
+        Gaussian-basis codes such as CRYSTAL (`tst`) and ORCA (`RMS-DP`).
+        """,
+    )
+
+    delta_density_matrix_max = Quantity(
+        shape=['*'],
+        type=float,
+        unit='dimensionless',
+        description="""
+        Maximum absolute change of the density-matrix elements `P_munu` between
+        consecutive SCF steps; the max-norm counterpart of
+        `delta_density_matrix_rms`. Reported by Gaussian-basis codes such as
+        CP2K, CRYSTAL (`PX`), and ORCA (`Max-DP`).
         """,
     )
 
@@ -280,55 +346,18 @@ def set_model_method_ref(self) -> ModelMethod | None:
 
     def _compute_energy_deltas(self, logger: 'BoundLogger'):
         """
-        Compute delta_energies_total from consecutive total_energies values.
-
-        Returns array of absolute energy differences between consecutive steps.
+        Compute `delta_energies_total` as the absolute change between consecutive
+        `scf_steps.energies_total` values, i.e. the per-SCF-iteration total-energy
+        series. This is the correct source for an SCF convergence measure; the
+        repeating `Outputs.total_energies` holds labeled generic total energies
+        (e.g. "DFT" vs "DFT + dispersion"), whose differences are not SCF deltas.
         """
-        try:
-            if self.total_energies is None or len(self.total_energies) < 2:
-                return None
-
-            # Extract energy values from each step
-            energy_values = [e.value for e in self.total_energies]
-
-            # Compute differences manually to preserve Pint units
-            deltas = []
-            for i in range(1, len(energy_values)):
-                delta = np.abs(energy_values[i] - energy_values[i - 1])
-                deltas.append(delta)
-
-            # Convert list to array-like structure
-            if len(deltas) > 0 and hasattr(deltas[0], 'magnitude'):
-                # Pint quantities - stack magnitudes and add units
-                magnitudes = [d.magnitude for d in deltas]
-                return np.array(magnitudes) * deltas[0].units
-            else:
-                return np.array(deltas)
-
-        except (AttributeError, IndexError, TypeError, ValueError) as e:
-            logger.debug(f'Could not compute delta_energies_total: {e}')
+        if self.scf_steps is None or self.scf_steps.energies_total is None:
             return None
-
-    def _compute_force_norms(self, logger: 'BoundLogger'):
-        """
-        Compute delta_force_abs from total_forces as force norms.
-
-        Returns array of force norms (L2 norm of 3D force vectors).
-        """
-        try:
-            if self.total_forces is None or len(self.total_forces) == 0:
-                return None
-
-            # Get force values (Pint Quantity with shape [n_atoms, 3])
-            force_values = self.total_forces[-1].value
-
-            # Compute force norms (preserves Pint units)
-            force_norms = ((force_values**2).sum(axis=1)) ** 0.5
-
-            return force_norms
-        except (AttributeError, IndexError, TypeError) as e:
-            logger.debug(f'Could not compute delta_force_abs: {e}')
+        energies = self.scf_steps.energies_total
+        if len(energies) < 2:
             return None
+        return np.abs(np.diff(energies.magnitude)) * energies.units
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
@@ -349,19 +378,20 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         except Exception as e:
             logger.debug(f'Could not set model_method_ref: {e}')
 
-        # Populate missing SCF delta quantities from available data
-        if self.scf_steps is not None:
-            # Define delta computations: (delta_field, compute_function)
-            delta_computations = [
-                ('delta_energies_total', self._compute_energy_deltas),
-                ('delta_force_abs', self._compute_force_norms),
-            ]
-
-            for delta_field, compute_func in delta_computations:
-                if getattr(self.scf_steps, delta_field, None) is None:
-                    computed_value = compute_func(logger)
-                    if computed_value is not None:
-                        setattr(self.scf_steps, delta_field, computed_value)
+        # Derive `delta_energies_total` from the per-SCF `energies_total` series when the
+        # parser did not provide it directly (#454).
+        #
+        # The other SCF residuals are deliberately NOT synthesized here. In particular
+        # `delta_force_abs` (the change of forces between successive SCF iterations) cannot
+        # be derived from the archive: `Outputs.total_forces` holds the final/labeled forces,
+        # not a per-SCF-iteration series, so norming them would fill an SCF-convergence
+        # field with the final force magnitudes -- a physically different quantity (#453).
+        # `delta_force_abs` and the density/potential residuals are therefore set only by
+        # parsers that genuinely report them per SCF step.
+        if self.scf_steps is not None and self.scf_steps.delta_energies_total is None:
+            deltas = self._compute_energy_deltas(logger)
+            if deltas is not None:
+                self.scf_steps.delta_energies_total = deltas
 
 
 class WorkflowOutputs(Outputs):