Real DFT as a level of theory, four more potentials, and two silent unit errors (#15)

* Real DFT as a level of theory, four more potentials, and two silent unit errors

matverse could reach six levels and every one of them was a surrogate. GPAW is
now registered as gpaw-pbe and gpaw-pbe-fast, with kind="dft" and
surrogate=False, so the one level that solves the Kohn-Sham equations is
distinguishable from the ones that reproduce something that did. m3gnet,
m3gnet-r2scan, tensornet and orb are added alongside.

Validated by running it, not by citing it. Through matverse, on real PBE:

  mv.calc.relax  -> Si lattice constant 5.479 A   (PBE literature 5.47)
  mv.prop.eos    -> Si bulk modulus     88.7 GPa  (PBE literature 88-89)
                    fit residual        0.002 meV/atom

Two things that would have been wrong silently.

matgl returns stress in GPa while ASE's contract is eV/A^3, so leaving
stress_unit alone makes every stress-derived quantity - elastic constants, bulk
modulus, pressure - exactly 160.2x too large with no error raised anywhere. The
factor is the GPa-to-eV/A^3 conversion, which is what makes it recognisable
once you look. Fixed, with a test that fails on the unfixed version.

The second took longer to find and is not GPAW's fault. I first gave the DFT
levels a k-point *density*, which is the sensible thing for comparing different
materials and the wrong thing for scanning one material's volume: the mesh
changes discretely as the cell grows, and each change puts a step in E(V) that
the fit reads as curvature. Silicon, same calculator, only the k-points varying:

  density 2.0 -> -879 GPa   (125 -> 64 k-points)
  density 2.5 ->  319 GPa   (216 -> 125)
  density 3.0 ->  125 GPa   (343 -> 216)
  density 4.0 -> 85.7 GPa   (729 -> 512)
  fixed 8x8x8 -> 88.7 GPa   (unchanged)

A negative bulk modulus is not a soft crystal. Raising the density only shrinks
the jump; it never removes it. The plane-wave cutoff was converged the whole
time - PW(400), PW(600) and PW(800) agree to 0.1 GPa - so the entire error was
the mesh. The built-in levels use a fixed mesh, and mv.prop.eos now says why.

obs['eos_residual'] catches it long before the modulus does: 9.5 meV/atom
against 0.002, a factor of nearly five thousand, and the advice already in that
function's notes - a residual far above a meV is a fit that should not be read -
was enough to reject it without knowing the cause.

On pymatviz I made a judgement rather than adding it wholesale. Most of it
duplicates mv.pl, which already has hull, parity, periodic_table, spectra,
scatter, bands and distribution; wrapping those again would be the mistake the
coverage map exists to prevent. The one genuinely absent capability is the
space-group distribution, so mv.pl.spacegroups does that directly, grouped and
coloured by crystal system rather than plotted as 230 flat categories.
mv.pp.symmetry gained spacegroup_number and spacegroup_symbol while I was
there - it was already building the analyser, and a symmetry report that gives
the crystal system and the point group but not the group itself is missing the
one label everybody cites.

960 tests pass on py3.12 with GPAW and matgl, 850 on py3.10 without them. The
level-metadata tests skip where a backend is absent rather than failing, which
they did not do at first. Contract-verified rate stays 165/165, all 16
notebooks execute, docs build.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016wh5NzNE9Wn49nuUm7zg2a

* Give GPAW and the torch-backed potentials a CI job of their own

All three test legs went red on the previous commit. GPAW has no wheel: it
compiles against libxc and BLAS, and the runner carries neither, so pip got as
far as "fatal error: xc.h: No such file or directory" before anything ran. It
built here only because this node happens to have both.

Installing it in the main matrix was the wrong call regardless of the headers.
GPAW compiles from source, and matgl and orb-models pull torch, which nothing
else in the matrix needs — three legs would each pay a couple of minutes to
cover one job's worth of tests. They now share a 'calculators' job, the same
shape the anharmonic job already has and for the same reason.

The job prints mv.calc.available() before running anything, so a level whose
backend failed to install is visible as a line in the log rather than as a
silently skipped test.

* Add build-essential for GPAW's C++ headers in the calculators job

GPAW 26 includes <algorithm> from sources the build compiles as C, so the
runner needs a C++ standard library it does not carry by default. The previous
attempt got past libxc and stopped at 'fatal error: algorithm: No such file or
directory'.

* Stop trying to build GPAW on the runner

Three attempts, same error. GPAW ships no wheel and its 26.x sources include
C++ headers from files the build compiles as C, so it stops at 'fatal error:
algorithm: No such file or directory' whether or not libxc, BLAS and
build-essential are installed. It builds here because this node's compiler puts
the C++ headers on the default include path; a stock ubuntu runner's does not.

So the calculators job installs the torch-backed potentials only, and keeps the
test that matters most there — the one that fails when matgl's stress is left
in GPa, a 160x error nothing else would catch.

The real-DFT numbers stay where they were measured, in the notes, and
TestRealDFT says why it skips rather than leaving a reader to assume CI covers
it.

---------

Co-authored-by: Zehua Zeng <steorra@sh04-04n05.int>
Co-authored-by: Claude Opus 5 <noreply@anthropic.com>

Author: 3 people

.github/workflows/test.yml

@@ -62,6 +62,50 @@ jobs:
         # failing assertion.
         run: pytest tests/test_contracts.py -q -s -k "rate or audit"
 
+  calculators:
+    # The torch-backed potentials get a job of their own: matgl and orb-models
+    # pull torch, which nothing else in this matrix needs, and making all three
+    # main legs carry it to cover one job's worth of tests is a bad trade.
+    #
+    # GPAW is deliberately NOT installed here. It ships no wheel, and its 26.x
+    # sources include C++ headers from files the build compiles as C, so on a
+    # stock runner it stops at "fatal error: algorithm: No such file or
+    # directory" whether or not libxc, BLAS and build-essential are present —
+    # three attempts, same error. It builds on machines whose compiler puts the
+    # C++ headers on the default include path, which is where the real-DFT
+    # numbers in the notes were measured: a silicon lattice constant of 5.479 A
+    # and a bulk modulus of 88.7 GPa. TestRealDFT skips here, and says so
+    # rather than pretending to cover it.
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
+      - name: Install
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e ".[dev,plot,potentials]"
+
+      - name: Show what is actually installed
+        run: |
+          python -c "
+          import matverse as mv
+          for name, meta in sorted(mv.calc.available().items()):
+              state = meta.get('unavailable', meta.get('method'))
+              print(f'{name:16s} {state}')
+          "
+
+      - name: Test the levels of theory
+        # MatglStress is the one that matters here: it fails on the version
+        # that leaves matgl's stress in GPa, which is a 160x error nothing
+        # else would catch.
+        run: |
+          pytest tests/test_pipeline.py -q \
+            -k "LevelsOfTheory or MatglStress" --durations=5
+
   anharmonic:
     # A job of its own, because hiphive reaches numba through trainstation and
     # numba caps numpy below 2.5. Installing it alongside the others would

matverse/calc.py

@@ -47,7 +47,9 @@
 _CALCULATORS: dict[str, tuple] = {}
 
 #: Levels matverse knows how to build, if their backend happens to be installed.
-BUILTIN_LEVELS = ("emt", "lj", "mace-mpa", "mace-omat", "sevennet", "chgnet")
+BUILTIN_LEVELS = ("emt", "lj", "mace-mpa", "mace-omat", "sevennet", "chgnet",
+                  "m3gnet", "m3gnet-r2scan", "tensornet", "orb",
+                  "gpaw-pbe", "gpaw-pbe-fast")
 
 
 @register_function(
@@ -122,6 +124,86 @@ def _builtin(name: str):
             "surrogate": True, "license": "BSD-3-Clause", "uncertainty": None,
             "note": "Superseded on Matbench Discovery by OMat24-trained models; "
                     "kept because a great deal of published work used it."})
+    if name in ("m3gnet", "m3gnet-r2scan", "tensornet"):
+        import matgl
+        from matgl.ext.ase import PESCalculator
+        checkpoint, functional = {
+            "m3gnet": ("M3GNet-PES-MatPES-PBE-2025.2", "PBE"),
+            "m3gnet-r2scan": ("M3GNet-PES-MatPES-r2SCAN-2025.2", "r2SCAN"),
+            "tensornet": ("TensorNet-PES-MatPES-PBE-2025.2", "PBE"),
+        }[name]
+        # stress_unit is not optional. matgl returns stress in GPa by
+        # default while ASE's contract is eV/A^3, so leaving it alone makes
+        # every stress-derived quantity - elastic constants, bulk modulus,
+        # pressure - exactly 160.2x too large, with no error anywhere.
+        return (lambda: PESCalculator(matgl.load_model(checkpoint),
+                                      stress_unit="eV/A3"),
+                {"kind": "mlip", "method": checkpoint,
+                 "reference": f"{functional} (MatPES)", "surrogate": True,
+                 "license": "BSD-3-Clause", "uncertainty": None,
+                 "note": "MatPES-trained, so it reproduces plain "
+                         f"{functional} rather than the PBE+U of the MP-"
+                         "trained generation. Which functional a surrogate "
+                         "was fitted to is not a detail: mixing a PBE+U "
+                         "surrogate with a PBE one is the same class of error "
+                         "as mixing PBE with HSE06, and it is why "
+                         "reference= is recorded on every level."})
+    if name == "orb":
+        from orb_models.forcefield import pretrained
+        from orb_models.forcefield.calculator import ORBCalculator
+        return (lambda: ORBCalculator(
+            pretrained.orb_v3_conservative_inf_omat(device="cpu"),
+            device="cpu"),
+            {"kind": "mlip", "method": "ORB v3 conservative", 
+             "reference": "PBE+U (OMat24)", "surrogate": True,
+             "license": "Apache-2.0", "uncertainty": None,
+             "note": "Conservative variant: forces are the energy gradient, "
+                     "which is what a phonon or a relaxation needs. The "
+                     "'direct' variants predict forces independently and are "
+                     "faster and not usable for either."})
+    if name in ("gpaw-pbe", "gpaw-pbe-fast"):
+        from gpaw import GPAW, PW
+        # Not a surrogate. This is the one level here that solves the
+        # Kohn-Sham equations rather than reproducing something that did, and
+        # kind/surrogate say so — every downstream record inherits it.
+        cutoff, mesh = ((500.0, (8, 8, 8)) if name == "gpaw-pbe"
+                        else (400.0, (6, 6, 6)))
+        return (lambda: GPAW(mode=PW(cutoff), xc="PBE", kpts=mesh,
+                             txt=None, symmetry={"point_group": False}),
+                {"kind": "dft", "method": f"GPAW PBE, PW({cutoff:g} eV)",
+                 "reference": "PBE", "surrogate": False,
+                 "license": "GPL-3.0", "uncertainty": None,
+                 "plane_wave_cutoff_eV": cutoff,
+                 "kpoint_mesh": list(mesh),
+                 "note": "Real plane-wave DFT, not a model of it. The cutoff "
+                         "and the k-point mesh are the two settings that "
+                         "decide whether a number is converged, so they are "
+                         "recorded on the level rather than left implicit.\n\n"
+                         "The mesh is **fixed** rather than set by a k-point "
+                         "density, and that is not a stylistic choice. A "
+                         "density-based mesh changes discretely as a cell "
+                         "changes size, which puts a step in E(V) and "
+                         "destroys anything fitted to it. Silicon, same "
+                         "calculator, only the k-points varying:\n\n"
+                         "  density 2.0 -> -879 GPa   (125 -> 64 k-points)\n"
+                         "  density 2.5 ->  319 GPa   (216 -> 125)\n"
+                         "  density 3.0 ->  125 GPa   (343 -> 216)\n"
+                         "  density 4.0 -> 85.7 GPa   (729 -> 512)\n"
+                         "  fixed 8x8x8 -> 88.7 GPa   (unchanged)\n\n"
+                         "against a PBE literature 88-89. A negative bulk "
+                         "modulus is not a soft crystal, it is a "
+                         "discontinuity. Raising the density only shrinks the "
+                         "relative size of the jump; it never removes it. The "
+                         "plane-wave cutoff, by contrast, was already "
+                         "converged - PW(400) and PW(600) differ by 0.1 GPa - "
+                         "so the whole error was the mesh.\n\n"
+                         "A fixed mesh suits small cells. Register your own "
+                         "level for anything large, where 8x8x8 is wasteful, "
+                         "and keep it fixed across any volume scan.\n\n"
+                         "point_group symmetry is off because matverse hands "
+                         "in displaced and strained cells whose symmetry is "
+                         "lower than the analyser infers from a rounded "
+                         "geometry."})
     raise KeyError(
         f"unknown level {name!r}. Runnable here: "
         f"{sorted(available(check_imports=False))}. Register your own with "

matverse/pl.py

@@ -17,6 +17,7 @@
 from __future__ import annotations
 
 import numpy as np
+import pandas as pd
 from anndata import AnnData
 
 from ._core import grid_of, structures
@@ -1162,3 +1163,99 @@ def distribution(md: AnnData, column: str, by: str | None = None,
         ax.spines[spine].set_visible(False)
     ax._matverse_dropped = dropped
     return ax
+
+
+@register_function(
+    aliases=["space group distribution", "symmetry distribution",
+             "spacegroup bar", "which space groups", "crystal system "
+             "distribution", "how symmetric is this set"],
+    category="pl",
+    description="Distribution of space groups across a dataset, grouped by "
+                "crystal system — what symmetry a generated or screened set "
+                "actually has.",
+    requires={"obs": ["{column}"]},
+    examples=["mv.pl.spacegroups(built)",
+              "mv.pl.spacegroups(md, column='spacegroup_number')",
+              "mv.pl.spacegroups(built, column='requested_space_group')"],
+    related=["mv.gen.from_symmetry", "mv.pp.describe", "mv.pl.distribution"],
+    notes="A generated set has a symmetry distribution, and it is rarely the "
+          "one that was asked for. mv.gen.from_symmetry records both the "
+          "requested group and the one the structure actually has, and "
+          "plotting the two side by side is how the difference becomes "
+          "visible rather than a column nobody reads.\\n\\n"
+          "Bars are grouped and coloured by crystal system rather than "
+          "plotted as 230 flat categories, because the number itself carries "
+          "no order a reader can use — 62 is not 'between' 61 and 63 in any "
+          "sense that matters, but Pnma being orthorhombic does.",
+)
+def spacegroups(md: AnnData, column: str = "spacegroup_number",
+                compare: str | None = None, top: int = 20, ax=None):
+    """Space-group distribution grouped by crystal system. Returns the axis."""
+    if column not in md.obs:
+        raise ValueError(
+            f"obs[{column!r}] absent; run mv.pp.symmetry(md) for "
+            f"'spacegroup_number', or point column= at "
+            f"mv.gen.from_symmetry's 'space_group' — this object has "
+            f"{sorted(md.obs.columns)[:8]}...")
+    if compare is not None and compare not in md.obs:
+        raise ValueError(f"obs[{compare!r}] absent")
+
+    numbers = pd.to_numeric(md.obs[column], errors="coerce").dropna()
+    if numbers.empty:
+        raise ValueError(f"obs[{column!r}] holds no usable space-group number")
+
+    counts = numbers.astype(int).value_counts().sort_values(ascending=False)
+    keep = counts.head(int(top))
+    order = sorted(keep.index)
+
+    ax = _axis(ax)
+    systems = [_crystal_system(n) for n in order]
+    palette = {"triclinic": "#8e44ad", "monoclinic": "#4c72b0",
+               "orthorhombic": "#2a9d8f", "tetragonal": "#e9c46a",
+               "trigonal": "#e76f51", "hexagonal": "#d35400",
+               "cubic": "#c1121f"}
+    positions = np.arange(len(order))
+    width = 0.4 if compare is not None else 0.7
+    ax.bar(positions - (width / 2 if compare is not None else 0),
+           [keep[n] for n in order], width=width,
+           color=[palette[s] for s in systems],
+           label=column if compare is not None else None)
+
+    if compare is not None:
+        other = pd.to_numeric(md.obs[compare], errors="coerce").dropna()
+        other = other.astype(int).value_counts()
+        ax.bar(positions + width / 2, [int(other.get(n, 0)) for n in order],
+               width=width, facecolor="none", edgecolor="#333333",
+               linewidth=0.9, label=compare)
+        ax.legend(frameon=False, fontsize=8)
+
+    ax.set_xticks(positions)
+    ax.set_xticklabels([str(n) for n in order], rotation=90, fontsize=7)
+    ax.set_xlabel("space group number")
+    ax.set_ylabel("materials")
+    for spine in ("top", "right"):
+        ax.spines[spine].set_visible(False)
+
+    seen = list(dict.fromkeys(systems))
+    handles = [_plt().Line2D([], [], color=palette[s], linewidth=6, label=s)
+               for s in seen]
+    ax.add_artist(ax.legend(handles=handles, frameon=False, fontsize=7,
+                            loc="upper right", title="crystal system",
+                            title_fontsize=7))
+    ax._matverse_n_groups = len(order)
+    ax._matverse_dropped = int(len(counts) - len(keep))
+    return ax
+
+
+#: Space-group number ranges, in the international convention.
+_CRYSTAL_SYSTEMS = ((2, "triclinic"), (15, "monoclinic"), (74, "orthorhombic"),
+                    (142, "tetragonal"), (167, "trigonal"), (194, "hexagonal"),
+                    (230, "cubic"))
+
+
+def _crystal_system(number: int) -> str:
+    """The crystal system a space-group number belongs to."""
+    for limit, name in _CRYSTAL_SYSTEMS:
+        if number <= limit:
+            return name
+    return "cubic"

matverse/pp.py

@@ -919,6 +919,7 @@ def prototype(md: AnnData, source: str = "input") -> None:
                 "which Wyckoff positions the atoms occupy.",
     requires={"structures": ["{source}"]},
     produces={"obs": ["crystal_system", "point_group",
+                      "spacegroup_number", "spacegroup_symbol",
                       "n_symmetry_operations", "n_wyckoff", "wyckoff",
                       "min_site_symmetry", "max_site_symmetry"]},
     examples=["mv.pp.symmetry(md)",
@@ -957,6 +958,8 @@ def symmetry(md: AnnData, source: str = "input", symprec: float = 0.01) -> None:
     systems = np.empty(md.n_obs, dtype=object)
     points = np.empty(md.n_obs, dtype=object)
     orders = np.full(md.n_obs, np.nan)
+    numbers = np.full(md.n_obs, np.nan)
+    group_symbols = np.empty(md.n_obs, dtype=object)
     wyckoff_count = np.full(md.n_obs, np.nan)
     wyckoff = np.empty(md.n_obs, dtype=object)
     lowest = np.full(md.n_obs, np.nan)
@@ -968,12 +971,19 @@ def symmetry(md: AnnData, source: str = "input", symprec: float = 0.01) -> None:
         systems[i] = ""
         points[i] = ""
         wyckoff[i] = ""
+        group_symbols[i] = ""
         try:
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
                 analyzer = SpacegroupAnalyzer(structure, symprec=symprec)
                 systems[i] = str(analyzer.get_crystal_system())
                 points[i] = str(analyzer.get_point_group_symbol())
+                # The analyser is already built, so the space group costs
+                # nothing extra - and a symmetry report that gives the crystal
+                # system and the point group but not the group itself is
+                # missing the one label everybody actually cites.
+                numbers[i] = float(analyzer.get_space_group_number())
+                group_symbols[i] = str(analyzer.get_space_group_symbol())
                 symbols = list(
                     analyzer.get_symmetrized_structure().wyckoff_symbols)
                 wyckoff[i] = ", ".join(str(w) for w in symbols)
@@ -993,6 +1003,8 @@ def symmetry(md: AnnData, source: str = "input", symprec: float = 0.01) -> None:
 
     md.obs["crystal_system"] = systems.astype(str)
     md.obs["point_group"] = points.astype(str)
+    md.obs["spacegroup_number"] = numbers
+    md.obs["spacegroup_symbol"] = group_symbols.astype(str)
     md.obs["n_symmetry_operations"] = orders
     md.obs["n_wyckoff"] = wyckoff_count
     md.obs["wyckoff"] = wyckoff.astype(str)

matverse/prop.py

@@ -229,7 +229,22 @@ def elastic(md: AnnData, level: str = "emt", source: str = "input",
           "absolute volume: materials of different size have no common volume "
           "axis, and the strain series they were computed on is common by "
           "construction. obs['eos_residual'] is the RMS misfit in eV/atom; a "
-          "value far above a meV is a fit that should not be read.",
+          "value far above a meV is a fit that should not be read.\n\n"
+          "**With a plane-wave DFT calculator, hold the k-point mesh fixed "
+          "across the scan.** A mesh set by a k-point *density* changes "
+          "discretely as the cell grows, and each change puts a step in E(V) "
+          "that the fit reads as curvature. On silicon with GPAW, varying "
+          "nothing but the k-points, the bulk modulus went -879, 319, 125 and "
+          "85.7 GPa at densities of 2.0, 2.5, 3.0 and 4.0, against 88.7 GPa "
+          "from the same calculator at a fixed 8x8x8 and a PBE literature "
+          "88-89. A negative bulk modulus is not a soft crystal; it is a "
+          "discontinuity. Raising the density only shrinks the jump.\n\n"
+          "obs['eos_residual'] catches it before the modulus does, and by a "
+          "wide margin: 9.5 meV/atom for the density-2.0 fit against 0.002 "
+          "for the fixed mesh, a factor of nearly five thousand. The existing "
+          "advice above — a residual far above a meV is a fit that should not "
+          "be read — was already enough to reject it. The built-in gpaw-pbe "
+          "levels use a fixed mesh so the situation does not arise.",
 )
 def eos(md: AnnData, level: str = "emt", source: str = "input",
         scales=None, model: str = "birch_murnaghan",

matverse_guide/docs/_scripts/nb_chemical_space.py

@@ -456,6 +456,27 @@ def mass_and_volume(structures):
 and everything else work on it unchanged. The funnel is now complete: elements →
 compositions → structures → energies.
 
+`mv.pl.spacegroups` shows the distribution, and putting the requested and the
+achieved group side by side is how the difference stops being a column nobody
+reads. Bars are grouped and coloured by crystal system rather than plotted as
+230 flat categories, because the number itself carries no order a reader can
+use — 62 is not "between" 61 and 63 in any sense that matters, but Pnma being
+orthorhombic does."""),
+
+    ("code", """\
+try:
+    ax = mv.pl.spacegroups(many, column="space_group",
+                           compare="requested_space_group")
+    ax.set_title("requested (hollow) vs achieved (filled)")
+except (ImportError, NameError, ValueError) as exc:
+    print("needs PyXtal and matplotlib:", exc)"""),
+
+    ("markdown", """\
+It also works on any dataset that has been through `mv.pp.symmetry`, which now
+records `spacegroup_number` and `spacegroup_symbol` alongside the crystal system
+and point group — the analyser was already being built, so the group itself cost
+nothing to report and is the one label everybody actually cites.
+
 ```{seealso}
 [Beyond one number](beyond_one_number.ipynb) covers the results that are not a
 single number per material: curves, per-atom values and measurements.