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

.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.