diff --git a/README.md b/README.md
index 2e80326..e0f4a8d 100644
--- a/README.md
+++ b/README.md
@@ -99,7 +99,7 @@ not single compositions.
 | `mv.gen` | scoring generated candidates, and enumerating substitutions |
 | `mv.model` | property prediction, with splits that do not leak |
 | `mv.opt` | design campaigns — what to compute next, and what came back |
-| `mv.pl` | plotting, including the periodic-table heatmap |
+| `mv.pl` | plotting — the periodic-table heatmap, hull, bands, phonon DOS and dispersion, Pourbaix map, NEB profile, MD RDF and MSD, Wulff shape, chemical-potential diagram |
 | `mv.utils` | units, checkpointing, cluster submission, object summaries |
 
 `mv.struct` is the v0.1 name for the structure half of `mv.pp`, kept as
@@ -267,7 +267,7 @@ Claims are deleted rather than repaired when they fail. The ones that did:
 | `utils.job_status requires uns['submissions']` | having submitted nothing is an answer, not an error |
 | `dft.status requires obs['dft_directory']` | scans the root directory, not the object |
 | `exp.attach requires uns['grids']` | a measured curve may be the first grid |
-| `surf.wulff produces uns['wulff']` | never written |
+| `surf.wulff produces uns['wulff']` | never written at the time. It is now — the polyhedron `mv.pl.wulff` draws lives there — and the claim is back as `bulk.uns['wulff']` |
 | `iface.build produces obs['nsites']` | never written |
 
 Two are informative beyond their own function. `mv.tl.cluster`'s two routes
diff --git a/matverse/md.py b/matverse/md.py
index 5c769f0..bc29957 100644
--- a/matverse/md.py
+++ b/matverse/md.py
@@ -137,6 +137,8 @@ def batched_available() -> dict:
     produces={"obs": ["md_energy_{level}", "md_temperature_{level}",
                       "msd_{level}", "diffusivity_{level}",
                       "md_volume_{level}"],
+              "obsm": ["md_temperature_trace_{level}", "md_msd_trace_{level}"],
+              "uns": ["grids"],
               "layers": ["diffusivity_{level}"],
               "structures": ["md_{level}"], "levels": ["{level}"]},
     prerequisites=["mv.calc.relax"],
@@ -150,7 +152,13 @@ def batched_available() -> dict:
           "the number that matters is the mobile species' diffusivity, not the "
           "average over everything in the cell, and averaging a lithium "
           "diffusivity with a framework one produces a number describing "
-          "nothing.",
+          "nothing.\n\n"
+          "Two traces are kept on the sampling grid: the temperature, so a "
+          "run that never equilibrated is visible rather than averaged away, "
+          "and the mean-squared displacement, which is the curve the "
+          "diffusivity is the slope of. mv.pl.rdf_msd draws the second one "
+          "with that slope on it; a diffusivity quoted without its MSD is a "
+          "number nobody can check.",
 )
 def run(md: AnnData, level: str = "emt", source: str = "input",
         temperature: float = 300.0, steps: int = 1000,
@@ -182,7 +190,7 @@ def run(md: AnnData, level: str = "emt", source: str = "input",
 
     energies, temperatures, msds, diffusivities = [], [], [], []
     volumes, finals, per_element, failed = [], [], [], 0
-    traces, times = [], None
+    traces, msd_traces, times = [], [], None
     elements = list(map(str, md.var_names))
     n_samples = len(range(0, steps, sample_every))
 
@@ -198,6 +206,7 @@ def run(md: AnnData, level: str = "emt", source: str = "input",
             volumes.append(np.nan); finals.append(structure)
             per_element.append(np.full(len(elements), np.nan))
             traces.append(np.full(n_samples, np.nan))
+            msd_traces.append(np.full(n_samples, np.nan))
             continue
         energies.append(result["energy"])
         temperatures.append(result["temperature"])
@@ -207,6 +216,7 @@ def run(md: AnnData, level: str = "emt", source: str = "input",
         finals.append(result["structure"])
         per_element.append(result["per_element"])
         traces.append(result["trace"])
+        msd_traces.append(result["msd_trace"])
         times = result["times"]
 
     md.obs[f"md_energy_{tag}"] = energies
@@ -218,6 +228,8 @@ def run(md: AnnData, level: str = "emt", source: str = "input",
         md.layers[f"diffusivity_{tag}"] = np.vstack(per_element)
     if times is not None and traces:
         _replace_trace(md, tag, np.vstack(traces), times)
+        _replace_trace(md, tag, np.vstack(msd_traces), times,
+                       quantity="md_msd_trace", unit="ps")
     deposit_structures(md, f"md_{tag}", finals)
     set_level(md, tag, **meta, source=source, ensemble=ensemble,
               temperature=temperature, steps=steps, timestep=timestep,
@@ -228,9 +240,10 @@ def run(md: AnnData, level: str = "emt", source: str = "input",
 
 
 def _replace_trace(md: AnnData, tag: str, block: np.ndarray,
-                   times: np.ndarray) -> None:
-    """Deposit the temperature trace, discarding any earlier one of a
-    different length.
+                   times: np.ndarray, quantity: str = "md_temperature_trace",
+                   unit: str = "K") -> None:
+    """Deposit a trace — temperature by default, MSD when asked — discarding
+    any earlier one of a different length.
 
     Grids exist so two levels of the same quantity can be compared, and
     ``deposit_grid`` refuses when a new grid disagrees with the stored one.
@@ -241,14 +254,13 @@ def _replace_trace(md: AnnData, tag: str, block: np.ndarray,
     so the stale one goes.
     """
     grids = md.uns.setdefault("grids", {})
-    stored = grids.get("md_temperature_trace")
+    stored = grids.get(quantity)
     if stored is not None and not np.array_equal(
             np.asarray(stored.get("values"), dtype=float), times):
-        for key in [k for k in md.obsm
-                    if k.startswith("md_temperature_trace_")]:
+        for key in [k for k in md.obsm if k.startswith(f"{quantity}_")]:
             del md.obsm[key]
-        del grids["md_temperature_trace"]
-    deposit_grid(md, "md_temperature_trace", tag, block, times, unit="K")
+        del grids[quantity]
+    deposit_grid(md, quantity, tag, block, times, unit=unit)
 
 
 def _check_thermostat(md: AnnData, tag: str, target: float,
@@ -393,6 +405,7 @@ def _integrate(structure, adaptor, calculator, temperature, steps, timestep,
             "temperature": kinetic_t / max(n, 1),
             "volume": volume / max(n, 1),
             "msd": float(squared.mean(axis=1)[-1]),
+            "msd_trace": squared.mean(axis=1),
             "diffusivity": overall,
             "per_element": by_element,
             "trace": np.asarray(trace, dtype=float),
diff --git a/matverse/pl.py b/matverse/pl.py
index 518851e..b81b667 100644
--- a/matverse/pl.py
+++ b/matverse/pl.py
@@ -1489,3 +1489,679 @@ def fermi_surface(md: AnnData, level: str = "dft", row=0,
         pass
     ax._matverse_n_sheets = len(sheets)
     return ax
+
+
+#: Okabe-Ito, the colour-blind-safe qualitative palette. Black is left out
+#: because reference lines - zero, the Fermi level, the water window - use it.
+OKABE_ITO = ("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2",
+             "#D55E00", "#CC79A7")
+
+
+def _row_labels(md: AnnData) -> list[str]:
+    """obs['name'] where a dataset carries one, obs_names otherwise."""
+    return [str(x) for x in md.obs.get("name", md.obs_names)]
+
+
+def _row_indices(md: AnnData, rows, limit: int = 5) -> list[int]:
+    if rows is None:
+        return list(range(min(md.n_obs, limit)))
+    return [int(r) for r in rows]
+
+
+def _resolve_row(md: AnnData, row) -> tuple[int, str, str]:
+    """``(index, obs_name, label)`` for a row given by name, obs_name or
+    position - the lookup order mv.pl.elastic settled on, because a matverse
+    dataset carries the formula in obs['name'] and integers in obs_names."""
+    labels, names = _row_labels(md), [str(x) for x in md.obs_names]
+    key = str(row)
+    if key in labels:
+        index = labels.index(key)
+    elif key in names:
+        index = names.index(key)
+    else:
+        try:
+            index = int(row)
+        except (TypeError, ValueError):
+            raise ValueError(f"no row {row!r}; this object has names "
+                             f"{labels[:8]} and an index of {names[:8]}") \
+                from None
+    return index, names[index], labels[index]
+
+
+@register_function(
+    aliases=["phonon plot", "plot phonon dos", "phonon density of states "
+             "plot", "plot phonons", "phonon dispersion plot", "vibrational "
+             "spectrum plot", "dispersion and dos"],
+    category="pl",
+    description="Draw the phonon density of states mv.prop.phonon stored, and "
+                "beside it the dispersion from mv.prop.dispersion when one is "
+                "passed, the two sharing the frequency axis.",
+    requires={"obsm": ["phonon_dos_{level}"], "uns": ["grids"]},
+    prerequisites=["mv.prop.phonon"],
+    examples=["mv.pl.phonon(md, level='emt')",
+              "ph = mv.prop.dispersion(md, level='emt')\n"
+              "mv.pl.phonon(md, level='emt', rows=[0], dispersion=ph)"],
+    related=["mv.prop.phonon", "mv.prop.dispersion", "mv.pl.bands",
+             "mv.pl.spectra"],
+    notes="mv.pl.spectra can overlay the DOS block like any other grid "
+          "quantity, and it would label the axis 'phonon_dos axis'. This "
+          "knows what it is drawing: frequency in THz, the DOS normalised to "
+          "unit area, the method and supercell the grid records in the "
+          "title, and the imaginary-mode count next to each material - "
+          "which is the one number the spectrum itself cannot show, because "
+          "mv.prop.phonon smears only the real modes onto the grid.\n\n"
+          "dispersion= takes the bands-axis object mv.prop.dispersion "
+          "returns. The two panels share the frequency axis, so the DOS is "
+          "drawn sideways - the layout every phonon code prints. The "
+          "dispersion panel is mv.pl.bands, and its abscissa is a fraction "
+          "along each material's own path rather than a wavevector; the "
+          "high-symmetry ticks are drawn only when one material is shown, "
+          "since two paths do not share them.\n\n"
+          "Returns the DOS axis. With a dispersion, the other panel is on "
+          "ax._matverse_dispersion_ax; pass ax=(left, right) to lay both "
+          "out yourself.",
+)
+def phonon(md: AnnData, level: str = "emt", rows=None, dispersion=None,
+           ax=None):
+    """Phonon DOS, with the dispersion beside it when given. Returns the axis."""
+    key = f"phonon_dos_{level}"
+    if key not in md.obsm:
+        raise ValueError(f"obsm[{key!r}] absent; run mv.prop.phonon(md, "
+                         f"level={level!r}) first")
+    grid = grid_of(md, "phonon_dos")
+    dos = np.asarray(md.obsm[key], dtype=float)
+    meta = md.uns.get("grids", {}).get("phonon_dos", {})
+    labels = _row_labels(md)
+    indices = _row_indices(md, rows)
+    imaginary = md.obs.get(f"n_imaginary_modes_{level}")
+
+    if dispersion is not None:
+        if ax is None:
+            figure = _plt().figure(figsize=(8.4, 4.2))
+            spec = figure.add_gridspec(1, 2, width_ratios=(3, 1), wspace=0.06)
+            left = figure.add_subplot(spec[0])
+            ax = figure.add_subplot(spec[1], sharey=left)
+        else:
+            left, ax = ax
+        known = list(dict.fromkeys(map(str, dispersion.obs["material"])))
+        names = [str(x) for x in md.obs_names]
+        wanted = [next(k for k in (names[i], labels[i]) if k in known)
+                  for i in indices if names[i] in known or labels[i] in known]
+        wanted = wanted or None
+        ticks = dispersion.uns.get("path_labels", {})
+        one = wanted[0] if wanted is not None and len(wanted) == 1 else None
+        bands(dispersion, materials=wanted, ax=left,
+              labels=ticks.get(one) if one is not None else None)
+        left.set_title(f"phonon dispersion ({level})", fontsize=10)
+        ax._matverse_dispersion_ax = left
+    else:
+        ax = _axis(ax)
+
+    counts = {}
+    for k, i in enumerate(indices):
+        label = labels[i]
+        n_im = int(imaginary.iloc[i]) if imaginary is not None else 0
+        counts[label] = n_im
+        if n_im > 0:
+            label += f" ({n_im} imaginary)"
+        elif n_im < 0:
+            label += " (failed)"
+        colour = OKABE_ITO[k % len(OKABE_ITO)]
+        if dispersion is not None:
+            ax.plot(dos[i], grid, color=colour, linewidth=1.0, label=label)
+            ax.fill_betweenx(grid, 0.0, dos[i], color=colour, alpha=0.15)
+        else:
+            ax.plot(grid, dos[i], color=colour, linewidth=1.0, label=label)
+            ax.fill_between(grid, 0.0, dos[i], color=colour, alpha=0.15)
+
+    parts = [str(meta.get("method") or "")]
+    if meta.get("supercell") is not None:
+        parts.append("x".join(str(int(s)) for s in meta["supercell"]))
+    detail = ", ".join(p for p in parts if p)
+    if dispersion is not None:
+        ax.set_xlabel("phonon DOS (1/THz)")
+        ax.tick_params(labelleft=False)
+        ax.set_title(f"DOS ({detail})" if detail else "DOS", fontsize=10)
+    else:
+        ax.set_xlabel("frequency (THz)")
+        ax.set_ylabel("phonon DOS (1/THz)")
+        ax.set_title(f"phonon DOS ({level}" + (f"; {detail}" if detail else "")
+                     + ")", fontsize=10)
+    ax.legend(frameon=False, fontsize=8)
+    for spine in ("top", "right"):
+        ax.spines[spine].set_visible(False)
+    ax._matverse_n_imaginary = counts
+    return ax
+
+
+@register_function(
+    aliases=["pourbaix diagram", "plot pourbaix", "aqueous stability plot",
+             "ph potential map", "corrosion diagram", "water stability plot",
+             "pourbaix map"],
+    category="pl",
+    description="Draw one material's aqueous decomposition energy over pH and "
+                "applied potential - the Pourbaix diagram as a stability "
+                "surface - with the water window and the point "
+                "mv.thermo.pourbaix evaluated marked on it.",
+    requires={"uns": ["pourbaix"]},
+    prerequisites=["mv.thermo.pourbaix"],
+    examples=["mv.pl.pourbaix(md)",
+              "mv.pl.pourbaix(md, row='LiFePO4', threshold=0.5)"],
+    related=["mv.thermo.pourbaix", "mv.thermo.hull", "mv.pl.hull"],
+    notes="A textbook Pourbaix diagram colours regions by which species "
+          "wins, which answers 'what does it become'. A screen asks a "
+          "different question - 'how far from stable is my candidate here' - "
+          "and that is a continuous quantity, the decomposition energy in "
+          "eV/atom, so it is drawn as a surface rather than as regions. The "
+          "white contour is the threshold below which a solid is commonly "
+          "taken to persist in water; 0.5 eV/atom is the value in use since "
+          "Singh et al. (2017), and it is an argument because it is a "
+          "convention rather than a law.\n\n"
+          "The dashed lines are the water stability window at 298 K - "
+          "oxygen evolution above, hydrogen evolution below. A candidate for "
+          "an aqueous electrode has to sit under the threshold between "
+          "them.\n\n"
+          "The map is read from uns['pourbaix']['maps'], which "
+          "mv.thermo.pourbaix stores per material because the diagram it "
+          "fetched to answer one point answers the whole plane at no extra "
+          "cost. A material it could not place has no map and is named in "
+          "uns['pourbaix']['errors'].",
+)
+def pourbaix(md: AnnData, row=0, threshold: float = 0.5,
+             cmap: str = "viridis", ax=None):
+    """Decomposition energy over pH and potential. Returns the axis."""
+    stored = md.uns.get("pourbaix")
+    if stored is None:
+        raise ValueError("uns['pourbaix'] absent; run mv.thermo.pourbaix(md, "
+                         "ph=..., potential=...) first")
+    maps = stored.get("maps") or {}
+    index, name, label = _resolve_row(md, row)
+    if name not in maps:
+        raise ValueError(
+            f"no Pourbaix map stored for {label!r}; mv.thermo.pourbaix keeps "
+            f"one per material it could place (here {sorted(maps)}), and a "
+            f"material it could not place is listed in "
+            f"uns['pourbaix']['errors']")
+    ph = np.asarray(stored["ph_grid"], dtype=float)
+    potential = np.asarray(stored["potential_grid"], dtype=float)
+    surface = np.asarray(maps[name], dtype=float)
+
+    ax = _axis(ax, figsize=(6.2, 4.4))
+    filled = ax.contourf(ph, potential, surface, levels=20, cmap=cmap)
+    ax.contour(ph, potential, surface, levels=[float(threshold)],
+               colors="white", linewidths=1.3)
+    # The water window at 298 K, in V vs SHE.
+    ax.plot(ph, 1.229 - 0.0591 * ph, color="#333333", linestyle="--",
+            linewidth=0.9)
+    ax.plot(ph, -0.0591 * ph, color="#333333", linestyle="--", linewidth=0.9)
+
+    at_ph, at_e = float(stored["ph"]), float(stored["potential"])
+    ax.scatter([at_ph], [at_e], marker="x", s=70, color="#D55E00",
+               linewidths=1.6, zorder=5)
+    column = "pourbaix_decomposition"
+    if column in md.obs:
+        value = float(md.obs[column].iloc[index])
+        ax.annotate(f"{value:.2f} eV/atom", (at_ph, at_e), xytext=(6, 6),
+                    textcoords="offset points", fontsize=8, color="#D55E00")
+
+    bar = ax.figure.colorbar(filled, ax=ax, pad=0.02)
+    bar.set_label("decomposition energy (eV/atom)")
+    ax.set_xlabel("pH")
+    ax.set_ylabel("E (V vs SHE)")
+    ax.set_title(f"{label} - aqueous stability; white contour at "
+                 f"{threshold:g} eV/atom", fontsize=10)
+    ax.set_xlim(float(ph.min()), float(ph.max()))
+    ax.set_ylim(float(potential.min()), float(potential.max()))
+    ax._matverse_threshold = float(threshold)
+    return ax
+
+
+@register_function(
+    aliases=["neb plot", "migration barrier plot", "energy profile",
+             "minimum energy path plot", "plot barrier", "reaction path plot",
+             "band profile"],
+    category="pl",
+    description="Draw the minimum-energy path mv.neb.barrier recorded - the "
+                "energy of each image against the path coordinate - with the "
+                "barrier marked and unconverged bands drawn dashed.",
+    requires={"obsm": ["neb_profile_{level}"], "uns": ["grids"]},
+    prerequisites=["mv.neb.barrier"],
+    examples=["mv.pl.neb(md, level='emt')",
+              "mv.pl.neb(md, level='emt', rows=[0, 2])"],
+    related=["mv.neb.barrier", "mv.neb.hop_endpoints", "mv.pl.spectra"],
+    notes="Images are drawn as points joined by lines rather than as a "
+          "curve, because a band of five images is five energies and a "
+          "smooth line through them would claim a resolution the "
+          "calculation does not have.\n\n"
+          "A band that did not converge is dashed and says so in the "
+          "legend. mv.neb.barrier records the number either way, and a "
+          "profile that has not converged still peaks somewhere; the dash "
+          "is the difference between a barrier and an artefact.\n\n"
+          "The shape is the check the number cannot make: one hump is a "
+          "single-hop mechanism, two humps an intermediate minimum and a "
+          "mechanism that was not modelled, and a profile that ends well "
+          "above zero a pair of endpoints that were not equivalent.",
+)
+def neb(md: AnnData, level: str = "emt", rows=None, ax=None):
+    """Energy profile along the band, barrier annotated. Returns the axis."""
+    key = f"neb_profile_{level}"
+    if key not in md.obsm:
+        raise ValueError(f"obsm[{key!r}] absent; run mv.neb.barrier(md, "
+                         f"initial=..., final=..., level={level!r}) first")
+    coordinate = grid_of(md, "neb_profile")
+    profiles = np.asarray(md.obsm[key], dtype=float)
+    converged = md.obs.get(f"neb_converged_{level}")
+    labels = _row_labels(md)
+
+    ax = _axis(ax)
+    barriers = {}
+    for k, i in enumerate(_row_indices(md, rows)):
+        row = profiles[i]
+        if not np.isfinite(row).any():
+            continue                     # the band failed; nothing to draw
+        ok = bool(converged.iloc[i]) if converged is not None else True
+        colour = OKABE_ITO[k % len(OKABE_ITO)]
+        ax.plot(coordinate, row, color=colour, linewidth=1.2, marker="o",
+                markersize=4.5, markeredgecolor="white",
+                linestyle="-" if ok else "--",
+                label=labels[i] + ("" if ok else " (not converged)"))
+        top = int(np.nanargmax(row))
+        barriers[labels[i]] = float(row[top])
+        ax.annotate(f"{row[top]:.2f} eV", (coordinate[top], row[top]),
+                    xytext=(0, 7), textcoords="offset points", ha="center",
+                    fontsize=8, color=colour)
+
+    ax.axhline(0.0, color="#333333", linewidth=0.8)
+    ax.set_xlim(float(coordinate.min()), float(coordinate.max()))
+    ax.set_xlabel("fractional path coordinate")
+    ax.set_ylabel("energy relative to the initial image (eV)")
+    ax.set_title(f"migration barrier ({level})", fontsize=10)
+    if barriers:
+        ax.legend(frameon=False, fontsize=8)
+    for spine in ("top", "right"):
+        ax.spines[spine].set_visible(False)
+    ax._matverse_barriers = barriers
+    return ax
+
+
+@register_function(
+    aliases=["rdf and msd", "md plot", "plot diffusion", "mean squared "
+             "displacement plot", "msd plot", "radial distribution plot",
+             "trajectory rdf plot", "diffusivity plot", "molecular dynamics "
+             "plot"],
+    category="pl",
+    description="Draw what a molecular dynamics run left behind: the "
+                "trajectory-averaged radial distribution function with its "
+                "running coordination number, and the mean-squared "
+                "displacement with the diffusivity's slope on it.",
+    requires={"obsm": ["rdf_md_{rdf_level}", "md_msd_trace_{level}"],
+              "uns": ["grids"]},
+    prerequisites=["mv.md.run", "mv.md.rdf"],
+    examples=["mv.pl.rdf_msd(md, level='emt')",
+              "mv.pl.rdf_msd(md, level='emt', which='msd')",
+              "mv.pl.rdf_msd(md, rdf_level='md', which='rdf')"],
+    related=["mv.md.run", "mv.md.rdf", "mv.md.conductivity", "mv.pl.spectra"],
+    notes="Two panels because they are the two things a diffusion claim "
+          "rests on. The RDF says whether the mobile species has a "
+          "structure at all - sharp shells are a solid, a smeared first "
+          "peak and nothing beyond it a liquid - and the running "
+          "coordination number (dotted, right axis) reads the shell "
+          "occupancy off it. The MSD is the curve the diffusivity is the "
+          "slope of, and the dashed line is that slope, 6Dt through the "
+          "second half of the run where mv.md.run fitted it. An MSD that "
+          "is still curving upward at the end has not reached the "
+          "diffusive regime, and the D on it is a fit to vibration.\n\n"
+          "The two come from different calls with different defaults: the "
+          "MSD from mv.md.run at level=, the RDF from mv.md.rdf at "
+          "rdf_level=, which is 'md' unless it was named. which='rdf' or "
+          "'msd' draws one panel alone into a single axis; the default "
+          "needs both. With both, the RDF axis is returned and the MSD "
+          "axis is on ax._matverse_msd_ax.",
+)
+def rdf_msd(md: AnnData, level: str = "emt", rdf_level: str = "md",
+            rows=None, which: str = "both", ax=None):
+    """RDF and MSD side by side. Returns the RDF axis (or the only axis)."""
+    if which not in ("both", "rdf", "msd"):
+        raise ValueError(f"which must be 'both', 'rdf' or 'msd', got {which!r}")
+    rdf_key, msd_key = f"rdf_md_{rdf_level}", f"md_msd_trace_{level}"
+    if which != "msd" and rdf_key not in md.obsm:
+        raise ValueError(f"obsm[{rdf_key!r}] absent; run mv.md.rdf(md, "
+                         f"trajectories, species=..., level={rdf_level!r}) "
+                         f"first, or pass which='msd'")
+    if which != "rdf" and msd_key not in md.obsm:
+        raise ValueError(f"obsm[{msd_key!r}] absent; run mv.md.run(md, "
+                         f"level={level!r}) first, or pass which='rdf'")
+    indices = _row_indices(md, rows)
+    labels = _row_labels(md)
+
+    if which == "both":
+        if ax is None:
+            _, (left, right) = _plt().subplots(1, 2, figsize=(10.0, 3.8),
+                                               layout="constrained")
+        else:
+            left, right = ax
+    else:
+        left = right = _axis(ax)
+
+    if which != "msd":
+        r = grid_of(md, "rdf_md")
+        g = np.asarray(md.obsm[rdf_key], dtype=float)
+        coordination_key = f"coordination_md_{rdf_level}"
+        twin = None
+        for k, i in enumerate(indices):
+            colour = OKABE_ITO[k % len(OKABE_ITO)]
+            left.plot(r, g[i], color=colour, linewidth=1.1, label=labels[i])
+            if coordination_key in md.obsm:
+                twin = twin if twin is not None else left.twinx()
+                twin.plot(r, np.asarray(md.obsm[coordination_key],
+                                        dtype=float)[i],
+                          color=colour, linewidth=0.9, linestyle=":")
+        shell = md.obs.get(f"first_shell_{rdf_level}")
+        if shell is not None:
+            for i in indices:
+                if np.isfinite(shell.iloc[i]):
+                    left.axvline(float(shell.iloc[i]), color="#999999",
+                                 linewidth=0.7, zorder=0)
+        if twin is not None:
+            twin.set_ylabel("running coordination number (dotted)")
+            twin.spines["top"].set_visible(False)
+        species = md.uns["grids"]["rdf_md"].get("species", "")
+        left.set_xlabel("r (Å)")
+        left.set_ylabel("g(r)")
+        left.set_title(f"{species} RDF over the trajectory ({rdf_level})"
+                       .strip(), fontsize=10)
+        left.legend(frameon=False, fontsize=8, loc="upper right")
+
+    if which != "rdf":
+        from .md import _A2_PER_PS_TO_CM2_PER_S
+
+        t = grid_of(md, "md_msd_trace")
+        msd = np.asarray(md.obsm[msd_key], dtype=float)
+        diffusivity = md.obs.get(f"diffusivity_{level}")
+        half = len(t) // 2
+        for k, i in enumerate(indices):
+            colour = OKABE_ITO[k % len(OKABE_ITO)]
+            label = labels[i]
+            d = float(diffusivity.iloc[i]) if diffusivity is not None \
+                else np.nan
+            if np.isfinite(d) and len(t) > half:
+                label += f"  D = {d:.2e} cm²/s"
+                slope = 6.0 * d / _A2_PER_PS_TO_CM2_PER_S       # Ų/ps
+                t0, y0 = t[half:].mean(), np.nanmean(msd[i][half:])
+                right.plot(t, y0 + slope * (t - t0), color=colour,
+                           linestyle="--", linewidth=0.9)
+            right.plot(t, msd[i], color=colour, linewidth=1.1, label=label)
+        right.set_xlabel("time (ps)")
+        right.set_ylabel("mean-squared displacement (Ų)")
+        right.set_title(f"MSD ({level}); dashed is 6Dt", fontsize=10)
+        right.legend(frameon=False, fontsize=8, loc="upper left")
+
+    for axis in {left, right}:
+        for spine in ("top", "right"):
+            axis.spines[spine].set_visible(False)
+    if which == "both":
+        left._matverse_msd_ax = right
+    return left
+
+
+@register_function(
+    aliases=["wulff plot", "plot wulff shape", "draw crystal shape",
+             "equilibrium shape plot", "nanoparticle shape plot", "show the "
+             "wulff construction", "particle shape"],
+    category="pl",
+    description="Draw the equilibrium crystal shape mv.surf.wulff built, in "
+                "three dimensions, one colour per Miller family with its share "
+                "of the surface.",
+    requires={"uns": ["wulff"]},
+    prerequisites=["mv.surf.wulff"],
+    examples=["mv.pl.wulff(md, level='emt')",
+              "mv.pl.wulff(md, level='emt', row='Cu', azimuth=20)"],
+    related=["mv.surf.wulff", "mv.surf.surface_energy", "mv.pl.fermi_surface"],
+    notes="Draws the polyhedron mv.surf.wulff kept in uns['wulff'] rather "
+          "than rebuilding it, so the plot shows exactly the construction "
+          "whose area fractions are on the facet rows.\n\n"
+          "There is no length axis on purpose. A Wulff shape has a shape "
+          "and no size: each face sits at a distance from the centre "
+          "proportional to its surface energy, so the coordinates are in "
+          "units of gamma and the particle is the same at ten nanometres "
+          "and ten microns. What the legend carries instead is the fraction "
+          "of the surface each family takes, which is the prediction an "
+          "electron micrograph can check.\n\n"
+          "A face is coloured by the Miller family of the plane that was "
+          "asked for, not by its symmetry copy - the eight (111) faces of a "
+          "cubic metal are one colour - and a family that mv.surf.wulff "
+          "found with zero area is absent here, because it is absent from "
+          "the particle.",
+)
+def wulff(md: AnnData, level: str = "emt", row=0, azimuth: float = 30.0,
+          elevation: float = 20.0, ax=None):
+    """The Wulff polyhedron in three dimensions. Returns the axis."""
+    stored = (md.uns.get("wulff") or {}).get(level)
+    if stored is None:
+        raise ValueError(f"uns['wulff'][{level!r}] absent; run "
+                         f"mv.surf.wulff(facets, bulk=md, level={level!r}) "
+                         f"first")
+    _, name, label = _resolve_row(md, row)
+    shape = stored.get(name)
+    if shape is None:
+        raise ValueError(
+            f"no Wulff shape was built for {label!r}; shapes exist for "
+            f"{sorted(stored)} - a material none of whose facets has a "
+            f"finite surface energy gets none")
+    if "vertices" not in shape:
+        raise ValueError("no geometry was stored for this shape; it was "
+                         "built by an older mv.surf.wulff - run it again")
+
+    from matplotlib.patches import Patch
+    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+
+    vertices = np.asarray(shape["vertices"], dtype=float)
+    face_index = np.asarray(shape["face_index"], dtype=int)
+    families = [str(f) for f in shape["face_miller"]]
+    kinds = list(dict.fromkeys(families))
+    colour = {k: OKABE_ITO[i % len(OKABE_ITO)] for i, k in enumerate(kinds)}
+
+    if ax is None:
+        ax = _plt().figure(figsize=(5.4, 5.0)).add_subplot(projection="3d")
+    polygons = [vertices[face_index == f] for f in range(len(families))]
+    _add_faces(ax, Poly3DCollection(
+        polygons, facecolors=[colour[f] for f in families],
+        edgecolors="#333333", linewidths=0.5, alpha=0.9))
+    extent = float(np.abs(vertices).max()) if vertices.size else 1.0
+    for setter in (ax.set_xlim, ax.set_ylim, ax.set_zlim):
+        setter(-extent, extent)
+
+    fractions = shape.get("area_fractions", {})
+    handles = [Patch(facecolor=colour[k], edgecolor="#333333",
+                     label=f"({' '.join(k.split('_'))})  "
+                           f"{100 * float(fractions.get(k, np.nan)):.0f}% "
+                           f"of the surface")
+               for k in kinds]
+    ax.legend(handles=handles, frameon=False, fontsize=8, loc="upper left")
+    ax.set_axis_off()
+    ax.set_title(f"{label} - Wulff shape ({level})\n"
+                 f"anisotropy {float(shape.get('anisotropy', np.nan)):.3f}",
+                 fontsize=10)
+    ax.view_init(elev=float(elevation), azim=float(azimuth))
+    try:
+        ax.set_box_aspect((1, 1, 1))
+    except Exception:                                      # pragma: no cover
+        pass
+    ax._matverse_n_faces = len(families)
+    return ax
+
+
+@register_function(
+    aliases=["chemical potential diagram plot", "plot chempot", "chempot "
+             "plot", "stability window plot", "phase stability region plot",
+             "growth conditions plot", "synthesis window plot"],
+    category="pl",
+    description="Draw the chemical potential diagram mv.thermo.chempot_diagram "
+                "stored - each phase's stability domain in the plane (binary) "
+                "or space (ternary) of elemental chemical potentials - or, with "
+                "kind='window', the per-phase ranges from "
+                "mv.thermo.chempot_limits.",
+    requires={"uns": ["chempot_diagram"]},
+    prerequisites=["mv.thermo.chempot_diagram"],
+    examples=["mv.pl.chempot(md)",
+              "mv.pl.chempot(md, limit=-3.0)",
+              "mv.pl.chempot(md, kind='window', element='O')"],
+    related=["mv.thermo.chempot_diagram", "mv.thermo.chempot_limits",
+             "mv.pl.hull"],
+    notes="In a system of n elements each phase's domain has dimension n-1: "
+          "a binary diagram is line segments in a plane, a ternary one is "
+          "polygons in a volume. Both are drawn from the vertices "
+          "mv.thermo.chempot_diagram stored; four or more elements have no "
+          "picture and the call says so rather than projecting one.\n\n"
+          "The elemental references have open domains that run to an "
+          "artificial floor - pymatgen's default_min_limit, -50 eV - which "
+          "is not a chemical potential anyone reaches. It is drawn at "
+          "limit=, one eV below the lowest physical vertex unless set, and "
+          "the title says where it was cut so a reader does not take the "
+          "edge of the picture for a boundary.\n\n"
+          "kind='window' is the other view of the same thing: for one "
+          "element, the range of its chemical potential over which each "
+          "phase stays on the hull, from mv.thermo.chempot_limits. On a hull "
+          "closed over one dataset the ranges are bounded by the dataset "
+          "rather than by chemistry, and the title carries that warning in "
+          "red.",
+)
+def chempot(md: AnnData, kind: str = "diagram", limit: float | None = None,
+            element: str | None = None, ax=None):
+    """Chemical potential domains (or windows). Returns the axis."""
+    if kind == "window":
+        return _chempot_window(md, element, ax)
+    if kind != "diagram":
+        raise ValueError(f"kind must be 'diagram' or 'window', got {kind!r}")
+    stored = md.uns.get("chempot_diagram")
+    if stored is None:
+        raise ValueError("uns['chempot_diagram'] absent; run "
+                         "mv.thermo.chempot_diagram(md, level=...) first")
+    elements = [str(e) for e in stored.get("elements", [])]
+    domains = {str(f): np.asarray(d["vertices"], dtype=float)
+               for f, d in stored["domains"].items() if len(d["vertices"])}
+    if len(elements) not in (2, 3) or not domains:
+        raise ValueError(
+            f"the diagram spans {len(elements)} elements ({elements}) with "
+            f"{len(domains)} domains; this draws a binary plane or a ternary "
+            f"volume - subset the object to a two- or three-element system")
+
+    every = np.vstack(list(domains.values()))
+    floor = float(every.min())                 # the artificial min limit
+    physical = every[every > floor + 1e-6]
+    display = float(limit) if limit is not None else \
+        (float(physical.min()) - 1.0 if physical.size else floor)
+    level = stored.get("level", "")
+
+    if len(elements) == 2:
+        ax = _axis(ax, figsize=(5.2, 5.0))
+        for k, (formula, points) in enumerate(domains.items()):
+            points = np.where(points <= floor + 1e-6, display, points)
+            order = np.lexsort((points[:, 1], points[:, 0]))
+            colour = OKABE_ITO[k % len(OKABE_ITO)]
+            ax.plot(points[order, 0], points[order, 1], color=colour,
+                    linewidth=3.0, solid_capstyle="round", label=formula)
+            ax.annotate(formula, points.mean(axis=0), fontsize=8,
+                        xytext=(4, 4), textcoords="offset points")
+        ax.set_xlim(display, 0.05 * abs(display))
+        ax.set_ylim(display, 0.05 * abs(display))
+        ax.set_xlabel(f"Δμ({elements[0]}) (eV)")
+        ax.set_ylabel(f"Δμ({elements[1]}) (eV)")
+        for spine in ("top", "right"):
+            ax.spines[spine].set_visible(False)
+    else:
+        from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+
+        if ax is None:
+            ax = _plt().figure(figsize=(5.8, 5.2)).add_subplot(projection="3d")
+        for k, (formula, points) in enumerate(domains.items()):
+            points = np.where(points <= floor + 1e-6, display, points)
+            if len(points) < 3:
+                continue
+            colour = OKABE_ITO[k % len(OKABE_ITO)]
+            _add_faces(ax, Poly3DCollection(
+                [_ring(points)], facecolor=colour, alpha=0.3,
+                edgecolor=colour, linewidths=0.8))
+            ax.text(*points.mean(axis=0), formula, fontsize=8)
+        for setter in (ax.set_xlim, ax.set_ylim, ax.set_zlim):
+            setter(display, 0.05 * abs(display))
+        ax.set_xlabel(f"Δμ({elements[0]}) (eV)")
+        ax.set_ylabel(f"Δμ({elements[1]}) (eV)")
+        ax.set_zlabel(f"Δμ({elements[2]}) (eV)")
+
+    ax.set_title(f"chemical potential diagram ({level}); open domains cut "
+                 f"at {display:.1f} eV", fontsize=9)
+    ax._matverse_n_domains = len(domains)
+    return ax
+
+
+def _add_faces(ax, collection) -> None:
+    """Add a 3D polygon collection with the axis limits left to the caller.
+
+    matplotlib >= 3.10 autoscales from the collection's vertex array, which
+    is NaN-padded when the faces have different vertex counts - as the faces
+    of a polyhedron do - and that path has produced infinite limits here.
+    Both callers set the limits themselves, so autoscaling is not needed;
+    older matplotlib has no such argument and never autoscaled.
+    """
+    try:
+        ax.add_collection3d(collection, autolim=False)
+    except TypeError:                                     # pragma: no cover
+        ax.add_collection3d(collection)
+
+
+def _ring(points: np.ndarray) -> np.ndarray:
+    """A planar polygon's vertices in perimeter order.
+
+    A domain is convex and flat, so ordering by angle about the centroid in
+    the plane's own basis gives the boundary; the basis comes from the two
+    leading singular vectors of the centred points.
+    """
+    centred = points - points.mean(axis=0)
+    _, _, basis = np.linalg.svd(centred, full_matrices=False)
+    flat = centred @ basis[:2].T
+    order = np.argsort(np.arctan2(flat[:, 1], flat[:, 0]))
+    return points[order]
+
+
+def _chempot_window(md: AnnData, element, ax):
+    """kind='window': per-phase chemical potential ranges for one element."""
+    stored = md.uns.get("chempot_limits")
+    if stored is None:
+        raise ValueError("uns['chempot_limits'] absent; run "
+                         "mv.thermo.chempot_limits(md, level=...) first")
+    limits = stored.get("limits", {})
+    seen = sorted({str(e) for ranges in limits.values() for e in ranges})
+    chosen = str(element) if element is not None else (seen[0] if seen else "")
+    rows = [(str(f), sorted(float(v) for v in r[chosen]))
+            for f, r in limits.items() if chosen in r]
+    if not rows:
+        raise ValueError(f"no window for element {chosen!r}; the limits "
+                         f"cover {seen}")
+    # Poorest in the element at the bottom, richest at the top, so the bars
+    # read as a ladder of growth conditions rather than in hull order.
+    rows.sort(key=lambda item: item[1])
+
+    ax = _axis(ax, figsize=(5.5, 0.4 * len(rows) + 1.2))
+    for k, (formula, (low, high)) in enumerate(rows):
+        ax.plot([low, high], [k, k], color=OKABE_ITO[k % len(OKABE_ITO)],
+                linewidth=6, solid_capstyle="butt")
+        ax.plot([low, high], [k, k], linestyle="none", marker="|",
+                color="#333333", markersize=11)
+    ax.set_yticks(range(len(rows)))
+    ax.set_yticklabels([f for f, _ in rows])
+    ax.set_xlabel(f"μ({chosen}) (eV)")
+    closed = bool(stored.get("closed_system", False))
+    ax.set_title(f"stability window ({stored.get('level', '')})"
+                 + (" - bounded by this dataset, not by chemistry"
+                    if closed else ""),
+                 fontsize=9, color="#c0392b" if closed else "black")
+    for spine in ("top", "right"):
+        ax.spines[spine].set_visible(False)
+    ax._matverse_n_windows = len(rows)
+    return ax
+
+
+__all__ += ["scatter", "bands", "distribution", "spacegroups", "elastic",
+            "fermi_surface", "phonon", "pourbaix", "neb", "rdf_msd", "wulff",
+            "chempot"]
diff --git a/matverse/surf.py b/matverse/surf.py
index af15daa..3e6d372 100644
--- a/matverse/surf.py
+++ b/matverse/surf.py
@@ -385,7 +385,8 @@ def surface_energy_chempot(facets: AnnData, bulk: AnnData, refs: AnnData,
     requires={"obs": ["surface_energy_{level}", "parent", "miller"]},
     produces={"facets.obs": ["wulff_area_fraction_{level}"],
               "bulk.obs": ["wulff_effective_radius_{level}",
-                           "wulff_shape_factor_{level}"]},
+                           "wulff_shape_factor_{level}"],
+              "bulk.uns": ["wulff"]},
     prerequisites=["mv.surf.surface_energy"],
     examples=["mv.surf.wulff(facets, bulk=md, level='emt')"],
     related=["mv.surf.surface_energy", "mv.surf.slabs"],
@@ -394,7 +395,13 @@ def surface_energy_chempot(facets: AnnData, bulk: AnnData, refs: AnnData,
           "an answer rather than a failure — that plane is not expressed.\n\n"
           "Deposits the per-facet area fractions onto the facet rows and the "
           "shape summary onto the bulk object, because the shape belongs to the "
-          "material and the fractions belong to its facets.",
+          "material and the fractions belong to its facets.\n\n"
+          "uns['wulff'][level][name] also keeps the polyhedron itself — the "
+          "vertices of every face on the shape, which face each belongs to, "
+          "and which Miller family each face is — so mv.pl.wulff can draw it "
+          "without rebuilding the construction. The vertex coordinates are in "
+          "the units of the surface energies (distance from the centre is "
+          "proportional to gamma), which is why the plot has no length axis.",
 )
 def wulff(facets: AnnData, bulk: AnnData, level: str = "emt",
           symprec: float = 0.1) -> None:
@@ -445,6 +452,9 @@ def wulff(facets: AnnData, bulk: AnnData, level: str = "emt",
             "expressed": [_miller_label(k) for k, v in areas.items()
                           if v > 1e-6],
             "anisotropy": float(shape.anisotropy),
+            "area_fractions": {_miller_label(k): float(v)
+                               for k, v in areas.items()},
+            **_polyhedron(shape),
         }
 
     facets.obs[f"wulff_area_fraction_{level}"] = fractions
@@ -455,6 +465,34 @@ def wulff(facets: AnnData, bulk: AnnData, level: str = "emt",
     record(facets, "surf.wulff", level=level)
 
 
+def _polyhedron(shape) -> dict:
+    """The faces of a Wulff shape as three h5ad-writable arrays.
+
+    Faces have different vertex counts, so a list of polygons would be ragged
+    and unwritable. Instead every vertex goes into one ``(n, 3)`` array in
+    face order, ``face_index`` says which face each vertex belongs to, and
+    ``face_miller`` names the Miller family of each face — enough to redraw
+    the polyhedron and colour it by family.
+    """
+    vertices, face_index, face_miller = [], [], []
+    # shape.facets holds every symmetry copy of every plane asked for;
+    # shape.on_wulff is per *family*. A copy that reached the surface has the
+    # vertices of the simplices lying on it, and one that did not has none.
+    for facet in shape.facets:
+        if not facet.points or not shape.on_wulff[facet.m_ind_orig]:
+            continue
+        ring = shape.get_line_in_facet(facet)
+        if len(ring) < 3:
+            continue
+        face_miller.append(_miller_label(shape.miller_list[facet.m_ind_orig]))
+        for point in ring:
+            vertices.append([float(c) for c in point])
+            face_index.append(len(face_miller) - 1)
+    return {"vertices": np.asarray(vertices, dtype=float).reshape(-1, 3),
+            "face_index": np.asarray(face_index, dtype=int),
+            "face_miller": face_miller}
+
+
 @register_function(
     aliases=["adsorption sites", "adsorbate", "adsorption", "binding sites",
              "place adsorbate", "surface sites"],
diff --git a/matverse/thermo.py b/matverse/thermo.py
index 4cb1447..ef2e136 100644
--- a/matverse/thermo.py
+++ b/matverse/thermo.py
@@ -566,17 +566,32 @@ def defect_formation(defective: AnnData, host: AnnData, level: str = "emt",
                 "whether it survives in water.",
     requires={"structures": ["input"]},
     produces={"obs": ["pourbaix_decomposition"], "uns": ["pourbaix"]},
-    examples=["mv.thermo.pourbaix(md, ph=7.0, potential=0.0)"],
-    related=["mv.thermo.hull"],
+    examples=["mv.thermo.pourbaix(md, ph=7.0, potential=0.0)",
+              "mv.thermo.pourbaix(md, ph=7.0, potential=0.0, "
+              "ph_range=(0, 14), n_grid=100)"],
+    related=["mv.thermo.hull", "mv.pl.pourbaix"],
     notes="Needs mp-api and an MP_API_KEY: aqueous stability is measured "
           "against the ion energies Materials Project fits, and there is no "
           "way to compute it from a candidate set alone. A material on the "
           "solid-state hull can still dissolve, which is why this is a "
-          "separate question rather than a column of the same one.",
+          "separate question rather than a column of the same one.\n\n"
+          "obs['pourbaix_decomposition'] is the number at the one (pH, E) "
+          "asked for. uns['pourbaix']['maps'][name] is the same quantity "
+          "over a grid of pH and potential — the Pourbaix diagram as a "
+          "decomposition-energy surface rather than as coloured regions, "
+          "since a screen wants to know how far from stable a candidate is, "
+          "not only which phase wins. The diagram is fetched once per "
+          "material either way, so the map costs nothing extra to keep, and "
+          "mv.pl.pourbaix draws it.",
 )
 def pourbaix(md: AnnData, ph: float = 7.0, potential: float = 0.0,
-             api_key: str | None = None) -> None:
-    """Distance from aqueous stability, in eV/atom, at one pH and potential."""
+             api_key: str | None = None, ph_range=(-2.0, 16.0),
+             potential_range=(-2.0, 3.0), n_grid: int = 60) -> None:
+    """Distance from aqueous stability, in eV/atom, at one pH and potential.
+
+    ``ph_range``, ``potential_range`` and ``n_grid`` set the map stored in
+    ``uns['pourbaix']['maps']``; the point value is unaffected by them.
+    """
     import os
 
     try:
@@ -593,9 +608,13 @@ def pourbaix(md: AnnData, ph: float = 7.0, potential: float = 0.0,
                          "energies and cannot be computed from candidates alone")
 
     S = structures(md, "input")
-    distances, failures = [], []
+    ph_grid = np.linspace(float(ph_range[0]), float(ph_range[1]), int(n_grid))
+    e_grid = np.linspace(float(potential_range[0]),
+                         float(potential_range[1]), int(n_grid))
+    ph_mesh, e_mesh = np.meshgrid(ph_grid, e_grid)
+    distances, maps, failures = [], {}, []
     with MPRester(key) as mpr:
-        for structure in S:
+        for name, structure in zip(map(str, md.obs_names), S):
             elements = sorted({str(el)
                                for el in structure.composition.elements})
             try:
@@ -610,6 +629,10 @@ def pourbaix(md: AnnData, ph: float = 7.0, potential: float = 0.0,
                     raise KeyError("no matching Pourbaix entry")
                 distances.append(float(diagram.get_decomposition_energy(
                     entry, pH=ph, V=potential)))
+                # pymatgen vectorises over pH and V, so the whole map is one
+                # call on the diagram already in hand.
+                maps[name] = np.asarray(diagram.get_decomposition_energy(
+                    entry, pH=ph_mesh, V=e_mesh), dtype=float)
             except Exception as exc:
                 distances.append(np.nan)
                 failures.append(f"{structure.composition.reduced_formula}: "
@@ -617,7 +640,9 @@ def pourbaix(md: AnnData, ph: float = 7.0, potential: float = 0.0,
 
     md.obs["pourbaix_decomposition"] = distances
     md.uns["pourbaix"] = {"ph": float(ph), "potential": float(potential),
-                          "n_failed": len(failures), "errors": failures[:10]}
+                          "n_failed": len(failures), "errors": failures[:10],
+                          "ph_grid": ph_grid, "potential_grid": e_grid,
+                          "maps": maps, "unit": "eV/atom"}
     record(md, "thermo.pourbaix", ph=ph, potential=potential)
 
 
diff --git a/matverse_guide/docs/Release_notes.md b/matverse_guide/docs/Release_notes.md
index 6fda801..97b0c77 100644
--- a/matverse_guide/docs/Release_notes.md
+++ b/matverse_guide/docs/Release_notes.md
@@ -1,5 +1,37 @@
 # Release notes
 
+## Unreleased
+
+### Six plots the analysis namespaces produced data for and `mv.pl` could not draw
+
+- `mv.pl.phonon` — the DOS from `mv.prop.phonon`, in THz, with the dispersion
+  from `mv.prop.dispersion` beside it on a shared frequency axis when passed.
+- `mv.pl.pourbaix` — aqueous decomposition energy over pH and potential, the
+  water window and the evaluated point. `mv.thermo.pourbaix` now stores that
+  map per material in `uns['pourbaix']['maps']` (`ph_range=`,
+  `potential_range=`, `n_grid=`); it costs nothing beyond the diagram it
+  already fetched.
+- `mv.pl.neb` — the minimum-energy path from `mv.neb.barrier`, barrier
+  annotated, unconverged bands dashed.
+- `mv.pl.rdf_msd` — the trajectory RDF from `mv.md.rdf` with its running
+  coordination number, and the MSD with the diffusivity's slope on it.
+  `mv.md.run` now keeps the MSD trace it fitted D to, as
+  `obsm['md_msd_trace_<level>']` on the same grid as the temperature trace,
+  and claims both traces.
+- `mv.pl.wulff` — the equilibrium shape in three dimensions, one colour per
+  Miller family with its share of the surface. `mv.surf.wulff` now stores the
+  polyhedron (`vertices`, `face_index`, `face_miller`, `area_fractions`) in
+  `uns['wulff'][level][name]`, and the `bulk.uns['wulff']` claim the README's
+  deleted-claim table recorded as never written is back, because it is.
+- `mv.pl.chempot` — the domains from `mv.thermo.chempot_diagram` as segments
+  (binary) or polygons (ternary), with the artificial floor cut where the title
+  says; `kind='window'` draws the per-phase ranges from
+  `mv.thermo.chempot_limits`.
+
+All six are registered with `requires` contracts that the probe battery
+checks, and each is called in a tutorial.
+
+
 ## v0.1.72
 
 **First PyPI release since v0.1.14.** The fifty-seven development versions in
diff --git a/matverse_guide/docs/_scripts/nb_defects_and_diffusion.py b/matverse_guide/docs/_scripts/nb_defects_and_diffusion.py
index 7411c44..6bba73e 100644
--- a/matverse_guide/docs/_scripts/nb_defects_and_diffusion.py
+++ b/matverse_guide/docs/_scripts/nb_defects_and_diffusion.py
@@ -390,9 +390,7 @@
 number."""),
 
     ("code", """\
-ax = mv.pl.spectra(copper, "neb_profile", levels=("emt",), rows=[0])
-ax.set_title("the minimum energy path")
-ax.set_ylabel("energy relative to start (eV)")"""),
+ax = mv.pl.neb(copper, level="emt")"""),
 
     ("markdown", """\
 The path rises to a single saddle and comes back down, which is what a
diff --git a/matverse_guide/docs/_scripts/nb_dynamics.py b/matverse_guide/docs/_scripts/nb_dynamics.py
index 09584c8..b10b7ca 100644
--- a/matverse_guide/docs/_scripts/nb_dynamics.py
+++ b/matverse_guide/docs/_scripts/nb_dynamics.py
@@ -365,6 +365,12 @@ def torchsim_factory():
 the coordination number comes out at **11.7 against a true 12** — the shortfall
 being the Gaussian smearing spilling past the cutoff.
 
+`mv.pl.rdf_msd` puts the two halves of a diffusion claim side by side: the RDF,
+with the running coordination number dotted on the right axis, and the
+mean-squared displacement `mv.md.run` kept, with the diffusivity's slope drawn
+through the half of the run it was fitted on. The copper cell above has both
+once it is given a trajectory to average — here the same stand-in, jittered.
+
 ```{note}
 That integral is computed here from its definition,
 $n(r) = \\int 4\\pi r^2 \\rho\\, g(r)\\, dr$, rather than taken from pymatgen's
@@ -374,6 +380,23 @@ def torchsim_factory():
 `first_shell_coordination` had better be one.
 ```"""),
 
+    ("code", """\
+cu_cell = mv.structures(copper, "input")[0]
+cu_frames = np.tile(np.array(cu_cell.frac_coords), (40, 1, 1))
+cu_frames = cu_frames + np.random.default_rng(0).normal(0, 0.01, cu_frames.shape)
+
+try:
+    mv.md.rdf(copper, cu_frames, species="Cu", reference="Cu", r_max=6.0)
+    ax = mv.pl.rdf_msd(copper, level="emt")
+except ImportError as exc:
+    print(exc)                 # pymatgen-analysis-diffusion is an extra"""),
+
+    ("markdown", """\
+An MSD that is still curving at the end of the run has not reached the
+diffusive regime, and the D written in the legend is then a fit to vibration.
+For copper at 300 K that is the right reading: the slope is a few 10⁻⁷ cm²/s
+of thermal rattling, not transport."""),
+
     ("markdown", """\
 ## How much of the cell do they visit?
 
diff --git a/matverse_guide/docs/_scripts/nb_getting_started.py b/matverse_guide/docs/_scripts/nb_getting_started.py
index c26d8c9..99284e3 100644
--- a/matverse_guide/docs/_scripts/nb_getting_started.py
+++ b/matverse_guide/docs/_scripts/nb_getting_started.py
@@ -304,6 +304,16 @@ def mace_factory():
 metals.obs[["name", "bulk_modulus_emt", "debye_temperature_emt",
             "thermal_conductivity_emt", "dynamically_stable_emt"]].round(1)"""),
 
+    ("markdown", """\
+`dynamically_stable_emt` is one bit per material. The spectrum it was read off
+is in `obsm["phonon_dos_emt"]`, and `mv.pl.phonon` draws it in THz with the
+method and supercell the grid recorded, so a stored spectrum says how coarse it
+is. Pass `dispersion=` the object `mv.prop.dispersion` returns and the two
+share the frequency axis."""),
+
+    ("code", """\
+ax = mv.pl.phonon(metals, level="emt")"""),
+
     ("markdown", """\
 Worth comparing against measured values rather than accepting:
 
diff --git a/matverse_guide/docs/_scripts/nb_screening.py b/matverse_guide/docs/_scripts/nb_screening.py
index 8d2f9cc..ebc1c3a 100644
--- a/matverse_guide/docs/_scripts/nb_screening.py
+++ b/matverse_guide/docs/_scripts/nb_screening.py
@@ -336,6 +336,15 @@ def cell(symbols):
 alni.obs[["formula", "energy_demo", "chempot_stable_demo",
           "chempot_window_demo"]].round(3)"""),
 
+    ("code", """\
+ax = mv.pl.chempot(alni)"""),
+
+    ("markdown", """\
+In a binary each domain is a line segment: the phase is stable only where its
+formation energy is exactly met, and the segment runs between the conditions
+where a neighbour takes over. The elemental references run off to an artificial
+floor, which is cut where the title says rather than drawn to −50 eV."""),
+
     ("markdown", """\
 `chempot_window` is the extent of the region in chemical potential space where
 each phase wins. **AlNi has three times the window of Al₃Ni**, and that is the
@@ -366,6 +375,7 @@ def cell(symbols):
 try:
     mv.thermo.pourbaix(md, ph=7.0, potential=0.0)
     print(md.obs[["formula", "pourbaix_decomposition"]].round(3))
+    ax = mv.pl.pourbaix(md)         # the map it stored, for the first row
 except (ValueError, ImportError) as exc:
     print(f"{type(exc).__name__}: {exc}")"""),
 
diff --git a/matverse_guide/docs/_scripts/nb_surfaces_and_adsorption.py b/matverse_guide/docs/_scripts/nb_surfaces_and_adsorption.py
index 16e9a87..fece341 100644
--- a/matverse_guide/docs/_scripts/nb_surfaces_and_adsorption.py
+++ b/matverse_guide/docs/_scripts/nb_surfaces_and_adsorption.py
@@ -227,7 +227,20 @@
 orders of magnitude larger, and flakes."""),
 
     ("code", """\
-bulk.uns["wulff"]["emt"]["0"]"""),
+{k: v for k, v in bulk.uns["wulff"]["emt"]["0"].items()
+ if k not in ("vertices", "face_index", "face_miller")}"""),
+
+    ("markdown", """\
+The summary leaves out three arrays: the vertices of every face on the
+particle, which face each belongs to, and each face's Miller family.
+`mv.surf.wulff` keeps them so the shape can be drawn without rebuilding the
+construction, and `mv.pl.wulff` is that drawing — one colour per family, with
+the share of the surface each takes in the legend. There is no length axis
+because a Wulff shape has none: faces sit at distances proportional to their
+surface energies, and the particle looks the same at any size."""),
+
+    ("code", """\
+ax = mv.pl.wulff(bulk, level="emt")"""),
 
     ("code", """\
 fig, ax = plt.subplots(figsize=(6, 3.6))
diff --git a/pyproject.toml b/pyproject.toml
index a074b0a..ec05e31 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -87,7 +87,12 @@ mlip = ["mace-torch"]                                # subset of 'potentials'
 # is installable only on a newer interpreter. mv.md.batched_available() reports
 # which situation an installation is in, and the ASE path works regardless.
 batched = ["torch-sim-atomistic"]
-dev = ["pytest", "build", "twine"]
+# bibtexparser: pymatgen reads the citations of its elemental price table with
+# the v1 API (`bibtexparser.bparser`), which v2 (2026) removed; pymatgen does
+# not pin it, so an unconstrained install on a fresh interpreter picks 2.x and
+# `mv.prop.cost` / prototype matching fail at import. Pin it where the test
+# environment is assembled.
+dev = ["pytest", "build", "twine", "bibtexparser<2"]
 
 [project.urls]
 Homepage = "https://github.com/omicverse/matverse"
diff --git a/tests/_contract_cases.py b/tests/_contract_cases.py
index 0bdd5b8..c5d804f 100644
--- a/tests/_contract_cases.py
+++ b/tests/_contract_cases.py
@@ -266,6 +266,49 @@ def vibrating():
     return md
 
 
+def aqueous():
+    """A dataset carrying the map mv.thermo.pourbaix deposits, for
+    mv.pl.pourbaix. The producer needs Materials Project, so the map is
+    written by hand in the producer's own shape."""
+    md = mv.data.from_compositions(["Fe2O3"])
+    ph, potential = np.linspace(-2, 16, 19), np.linspace(-2, 3, 11)
+    P, E = np.meshgrid(ph, potential)
+    md.obs["pourbaix_decomposition"] = [0.3]
+    md.uns["pourbaix"] = {
+        "ph": 7.0, "potential": 0.0, "n_failed": 0, "errors": [],
+        "ph_grid": ph, "potential_grid": potential, "unit": "eV/atom",
+        "maps": {"0": 0.6 * np.abs(E - (0.8 - 0.059 * P))}}
+    return md
+
+
+def barred():
+    """A band already run, for mv.pl.neb."""
+    md = hopped()
+    mv.neb.barrier(md, "hop_initial", "hop_final", level="emt", n_images=3,
+                   steps=5)
+    return md
+
+
+def diffusing():
+    """An MD run and a trajectory RDF on one cell, for mv.pl.rdf_msd. The RDF
+    needs pymatgen-analysis-diffusion; without it the probe is undecided."""
+    md = mv.datasets.metals(["Cu"], supercell=(2, 2, 2))
+    mv.pp.describe(md)
+    mv.md.run(md, level="emt", steps=30, equilibration=10, sample_every=5)
+    cell = mv.structures(md, "input")[0]
+    frames = np.tile(np.array(cell.frac_coords), (10, 1, 1))
+    frames = frames + np.random.default_rng(0).normal(0, 0.01, frames.shape)
+    mv.md.rdf(md, frames, species="Cu", reference="Cu", r_max=6.0)
+    return md
+
+
+def chempot_mapped():
+    """Domains already computed, for mv.pl.chempot."""
+    md = with_alni()
+    mv.thermo.chempot_diagram(md, level="emt")
+    return md
+
+
 def adsorbing():
     """Synthetic binding energies across a set of surfaces."""
     n = 6
@@ -655,6 +698,12 @@ def wulffable():
         mv.surf.surface_energy(out, bulk, level="emt")
         return out
 
+    def wulffed():
+        """The bulk after mv.surf.wulff, for mv.pl.wulff."""
+        out = bulk.copy()
+        mv.surf.wulff(wulffable(), out, level="emt")
+        return out
+
     # The off-stoichiometry route needs slabs that are actually off it, and an
     # elemental structure has none: every cut of Cu is Cu. Ordered Cu3Au cut
     # with symmetrize=True is the smallest case that produces a surface excess,
@@ -1151,6 +1200,12 @@ def measured_alloy():
         (mv.pl.spectra, patterned, ("xrd",), {}),
         (mv.pl.provenance, described, (), {}),
         (mv.pl.rank_elements_groups, ranked, (), {}),
+        (mv.pl.phonon, vibrating, (), {"level": "emt"}),
+        (mv.pl.pourbaix, aqueous, (), {}),
+        (mv.pl.neb, barred, (), {"level": "emt"}),
+        (mv.pl.rdf_msd, diffusing, (), {"level": "emt"}),
+        (mv.pl.wulff, wulffed, (), {"level": "emt"}),
+        (mv.pl.chempot, chempot_mapped, (), {}),
     ]
 
 
diff --git a/tests/test_pl_chempot.py b/tests/test_pl_chempot.py
new file mode 100644
index 0000000..b39a812
--- /dev/null
+++ b/tests/test_pl_chempot.py
@@ -0,0 +1,106 @@
+"""mv.pl.chempot: the domains mv.thermo.chempot_diagram stored, and the
+windows mv.thermo.chempot_limits stored."""
+
+from __future__ import annotations
+
+import pytest
+
+import matverse as mv
+
+
+@pytest.fixture(autouse=True)
+def _agg():
+    import matplotlib
+    matplotlib.use("Agg")
+
+
+def _fcc(symbol, a):
+    from pymatgen.core import Lattice, Structure
+    return Structure(Lattice.cubic(a), [symbol] * 4,
+                     [[0, 0, 0], [0, .5, .5], [.5, 0, .5], [.5, .5, 0]])
+
+
+def _b2(a_sym, b_sym, a):
+    from pymatgen.core import Lattice, Structure
+    return Structure(Lattice.cubic(a), [a_sym, b_sym],
+                     [[0, 0, 0], [.5, .5, .5]])
+
+
+def _l12(host, guest, a):
+    from pymatgen.core import Lattice, Structure
+    return Structure(Lattice.cubic(a), [guest, host, host, host],
+                     [[0, 0, 0], [0, .5, .5], [.5, 0, .5], [.5, .5, 0]])
+
+
+@pytest.fixture(scope="module")
+def binary():
+    """Al-Ni with formation energies supplied: Al3Ni at -1.8, AlNi at -1.4."""
+    md = mv.data.from_structures([_l12("Al", "Ni", 3.78), _b2("Al", "Ni", 2.89),
+                                  _fcc("Al", 4.05), _fcc("Ni", 3.52)])
+    mv.pp.describe(md)
+    md.obs["energy_demo"] = [-1.8, -1.4, 0.0, 0.0]
+    mv.thermo.chempot_diagram(md, level="demo")
+    mv.thermo.chempot_limits(md, level="demo")
+    return md
+
+
+@pytest.fixture(scope="module")
+def ternary():
+    md = mv.data.from_structures([_fcc("Al", 4.05), _fcc("Cu", 3.61),
+                                  _fcc("Ni", 3.52), _b2("Al", "Ni", 2.89)])
+    mv.calc.energy(md, level="emt")
+    mv.thermo.chempot_diagram(md, level="emt")
+    return md
+
+
+class TestDiagram:
+    def test_a_binary_is_segments_in_the_plane(self, binary):
+        ax = mv.pl.chempot(binary)
+        assert ax._matverse_n_domains == 4
+        assert len(ax.lines) == 4
+        assert ax.get_xlabel() == "Δμ(Al) (eV)"
+        assert ax.get_ylabel() == "Δμ(Ni) (eV)"
+        assert "demo" in ax.get_title() and "cut at" in ax.get_title()
+        # The artificial floor at -50 eV is not drawn; the axis stops one eV
+        # below the lowest physical vertex (-1.8).
+        assert ax.get_xlim()[0] == pytest.approx(-2.8)
+        for line in ax.lines:
+            assert (line.get_xdata() >= -2.8 - 1e-9).all()
+
+    def test_the_floor_is_an_argument(self, binary):
+        ax = mv.pl.chempot(binary, limit=-4.0)
+        assert ax.get_ylim()[0] == pytest.approx(-4.0)
+
+    def test_a_ternary_is_polygons_in_a_volume(self, ternary):
+        ax = mv.pl.chempot(ternary)
+        assert ax.name == "3d"
+        assert ax._matverse_n_domains >= 3
+        assert len(ax.collections) == ax._matverse_n_domains
+        assert ax.get_zlabel().startswith("Δμ(")
+
+    def test_says_what_is_missing(self):
+        md = mv.data.from_compositions(["AlNi"])
+        with pytest.raises(ValueError, match="mv.thermo.chempot_diagram"):
+            mv.pl.chempot(md)
+        with pytest.raises(ValueError, match="kind must be"):
+            mv.pl.chempot(md, kind="triangle")
+
+
+class TestWindow:
+    def test_one_bar_per_phase_with_a_window(self, binary):
+        ax = mv.pl.chempot(binary, kind="window", element="Ni")
+        # The elemental references have no window; the two compounds do.
+        assert ax._matverse_n_windows == 2
+        assert [t.get_text() for t in ax.get_yticklabels()] == ["Al3Ni", "AlNi"]
+        assert ax.get_xlabel() == "μ(Ni) (eV)"
+        assert "bounded by this dataset" in ax.get_title()
+
+    def test_the_element_defaults_to_the_first_seen(self, binary):
+        assert mv.pl.chempot(binary, kind="window").get_xlabel() == "μ(Al) (eV)"
+
+    def test_says_what_is_missing(self, binary):
+        with pytest.raises(ValueError, match="no window for element 'O'"):
+            mv.pl.chempot(binary, kind="window", element="O")
+        md = mv.data.from_compositions(["AlNi"])
+        with pytest.raises(ValueError, match="mv.thermo.chempot_limits"):
+            mv.pl.chempot(md, kind="window")
diff --git a/tests/test_pl_neb.py b/tests/test_pl_neb.py
new file mode 100644
index 0000000..cde0e19
--- /dev/null
+++ b/tests/test_pl_neb.py
@@ -0,0 +1,58 @@
+"""mv.pl.neb: the profile mv.neb.barrier recorded, with the barrier on it."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import matverse as mv
+
+
+@pytest.fixture(autouse=True)
+def _agg():
+    import matplotlib
+    matplotlib.use("Agg")
+
+
+@pytest.fixture(scope="module")
+def barred():
+    md = mv.datasets.metals(["Cu"])
+    mv.pp.describe(md)
+    mv.neb.hop_endpoints(md, "Cu", supercell=(1, 1, 1), key_added="hop")
+    mv.neb.barrier(md, "hop_initial", "hop_final", level="emt", n_images=5,
+                   steps=5)
+    return md
+
+
+class TestNeb:
+    def test_draws_the_images_and_annotates_the_barrier(self, barred):
+        ax = mv.pl.neb(barred, level="emt")
+        (line,) = [l for l in ax.lines if len(l.get_xdata()) == 5]
+        assert line.get_marker() == "o", "images are points, not a curve"
+        np.testing.assert_allclose(line.get_xdata(), np.linspace(0, 1, 5))
+        assert "eV" in ax.get_ylabel()
+        assert "path coordinate" in ax.get_xlabel()
+        assert "emt" in ax.get_title()
+        barrier = float(barred.obs["barrier_emt"].iloc[0])
+        assert ax._matverse_barriers == {"Cu": pytest.approx(barrier)}
+        assert any(f"{barrier:.2f} eV" == t.get_text() for t in ax.texts)
+
+    def test_an_unconverged_band_is_dashed_and_says_so(self, barred):
+        md = barred.copy()
+        md.obs["neb_converged_emt"] = [False]
+        ax = mv.pl.neb(md, level="emt")
+        (line,) = [l for l in ax.lines if len(l.get_xdata()) == 5]
+        assert line.get_linestyle() == "--"
+        assert "not converged" in line.get_label()
+
+    def test_a_failed_band_is_skipped_rather_than_drawn(self, barred):
+        md = barred.copy()
+        md.obsm["neb_profile_emt"] = np.full_like(md.obsm["neb_profile_emt"],
+                                                  np.nan)
+        ax = mv.pl.neb(md, level="emt")
+        assert ax._matverse_barriers == {}
+
+    def test_says_what_is_missing(self):
+        md = mv.datasets.metals(["Cu"])
+        with pytest.raises(ValueError, match="mv.neb.barrier"):
+            mv.pl.neb(md, level="emt")
diff --git a/tests/test_pl_phonon.py b/tests/test_pl_phonon.py
new file mode 100644
index 0000000..3486e58
--- /dev/null
+++ b/tests/test_pl_phonon.py
@@ -0,0 +1,77 @@
+"""mv.pl.phonon: the DOS mv.prop.phonon stored, and the dispersion beside it."""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+from anndata import AnnData
+
+import matverse as mv
+
+
+@pytest.fixture(autouse=True)
+def _agg():
+    import matplotlib
+    matplotlib.use("Agg")
+
+
+@pytest.fixture(scope="module")
+def vibrating():
+    md = mv.datasets.metals(["Cu", "Al"])
+    mv.pp.describe(md)
+    mv.prop.phonon(md, level="emt", supercell=(1, 1, 1))
+    return md
+
+
+def _dispersion_like(names, n_points=20, n_branches=3):
+    """The layout mv.prop.dispersion returns, built without phonopy."""
+    rows, material = [], []
+    for name in names:
+        for b in range(n_branches):
+            rows.append(np.abs(np.sin(np.linspace(0, np.pi, n_points))) * (b + 1))
+            material.append(name)
+    ph = AnnData(
+        X=np.vstack(rows),
+        obs=pd.DataFrame({"material": pd.Categorical(material)},
+                         index=[f"{m}:{i}" for i, m in enumerate(material)]),
+        var=pd.DataFrame({"path_fraction": np.linspace(0, 1, n_points)},
+                         index=[f"q{i}" for i in range(n_points)]))
+    ph.uns["y_label"] = "frequency (THz)"
+    ph.uns["path_labels"] = {name: {0.0: "Γ", 0.5: "X", 1.0: "Γ"} for name in names}
+    return ph
+
+
+class TestPhonon:
+    def test_draws_one_curve_per_row_with_units(self, vibrating):
+        ax = mv.pl.phonon(vibrating, level="emt")
+        assert len(ax.lines) == 2
+        assert "THz" in ax.get_xlabel()
+        assert "1/THz" in ax.get_ylabel()
+        assert "emt" in ax.get_title()
+        assert "commensurate" in ax.get_title(), "the grid records the method"
+        assert set(ax._matverse_n_imaginary) == {"Cu", "Al"}
+
+    def test_the_dispersion_shares_the_frequency_axis(self, vibrating):
+        ph = _dispersion_like([str(n) for n in vibrating.obs_names])
+        ax = mv.pl.phonon(vibrating, level="emt", rows=[0], dispersion=ph)
+        left = ax._matverse_dispersion_ax
+        assert left.get_shared_y_axes().joined(left, ax)
+        assert left._matverse_n_bands == 3, "only the requested row's branches"
+        assert [t.get_text() for t in left.get_xticklabels()] == ["Γ", "X", "Γ"]
+        assert "1/THz" in ax.get_xlabel()
+        assert len(ax.lines) == 1
+
+    def test_imaginary_modes_are_named_in_the_legend(self, vibrating):
+        md = vibrating.copy()
+        md.obs["n_imaginary_modes_emt"] = [2, 0]
+        ax = mv.pl.phonon(md, level="emt")
+        texts = [t.get_text() for t in ax.get_legend().get_texts()]
+        assert texts[0].endswith("(2 imaginary)") and texts[1] == "Al"
+
+    def test_says_what_is_missing(self, vibrating):
+        bare = mv.datasets.metals(["Cu"])
+        with pytest.raises(ValueError, match="mv.prop.phonon"):
+            mv.pl.phonon(bare, level="emt")
+        with pytest.raises(ValueError, match="mv.prop.phonon"):
+            mv.pl.phonon(vibrating, level="not-run")
diff --git a/tests/test_pl_pourbaix.py b/tests/test_pl_pourbaix.py
new file mode 100644
index 0000000..7b1f618
--- /dev/null
+++ b/tests/test_pl_pourbaix.py
@@ -0,0 +1,85 @@
+"""mv.pl.pourbaix draws the map mv.thermo.pourbaix stores.
+
+The producer needs Materials Project's aqueous entries and a key, so the map
+here is synthetic and shaped exactly as mv.thermo.pourbaix deposits it. The
+plot is what is under test; the deposit's shape is pinned by the probe of the
+producer's own claim.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import matverse as mv
+
+
+@pytest.fixture(autouse=True)
+def _agg():
+    import matplotlib
+    matplotlib.use("Agg")
+
+
+@pytest.fixture
+def aqueous():
+    md = mv.data.from_compositions(["Fe2O3", "TiO2"])
+    md.obs["name"] = ["Fe2O3", "TiO2"]
+    ph = np.linspace(-2.0, 16.0, 37)
+    potential = np.linspace(-2.0, 3.0, 26)
+    P, E = np.meshgrid(ph, potential)
+    md.obs["pourbaix_decomposition"] = [0.35, np.nan]
+    md.uns["pourbaix"] = {
+        "ph": 7.0, "potential": 0.0, "n_failed": 1,
+        "errors": ["TiO2: KeyError: no matching Pourbaix entry"],
+        "ph_grid": ph, "potential_grid": potential, "unit": "eV/atom",
+        # A band of stability along a Nernst line, as an oxide would show.
+        "maps": {"0": 0.6 * np.abs(E - (0.8 - 0.059 * P))},
+    }
+    return md
+
+
+class TestPourbaix:
+    def test_draws_the_surface_the_window_and_the_point(self, aqueous):
+        ax = mv.pl.pourbaix(aqueous)
+        assert ax.collections, "the filled contours"
+        assert len(ax.lines) == 2, "the two water-window lines"
+        assert ax.get_xlabel() == "pH"
+        assert "V vs SHE" in ax.get_ylabel()
+        labels = [t.get_text() for t in ax.texts]
+        assert any("0.35 eV/atom" in t for t in labels), \
+            "the value at the evaluated point is written next to it"
+        assert "Fe2O3" in ax.get_title()
+        colorbar = ax.figure.axes[-1]
+        assert "eV/atom" in colorbar.get_ylabel()
+
+    def test_the_threshold_is_an_argument(self, aqueous):
+        ax = mv.pl.pourbaix(aqueous, threshold=0.2)
+        assert ax._matverse_threshold == 0.2
+        assert "0.2" in ax.get_title()
+
+    def test_a_row_without_a_map_is_named(self, aqueous):
+        with pytest.raises(ValueError, match="no Pourbaix map stored for 'TiO2'"):
+            mv.pl.pourbaix(aqueous, row="TiO2")
+        with pytest.raises(ValueError, match="no Pourbaix map stored for 'TiO2'"):
+            mv.pl.pourbaix(aqueous, row=1)
+
+    def test_says_what_is_missing(self):
+        md = mv.data.from_compositions(["Fe2O3"])
+        with pytest.raises(ValueError, match="mv.thermo.pourbaix"):
+            mv.pl.pourbaix(md)
+
+
+class TestProducerSignature:
+    def test_the_producer_takes_the_map_grid(self):
+        """The map is deposited by mv.thermo.pourbaix; without a key it must
+        still refuse before touching the network, with the new arguments
+        accepted."""
+        import inspect
+        params = inspect.signature(mv.thermo.pourbaix).parameters
+        assert {"ph_range", "potential_range", "n_grid"} <= set(params)
+        md = mv.data.from_compositions(["Fe2O3"])
+        try:
+            with pytest.raises(ValueError, match="MP_API_KEY"):
+                mv.thermo.pourbaix(md, api_key=None, n_grid=5)
+        except ImportError:
+            pytest.skip("mp-api is not installed")
diff --git a/tests/test_pl_rdf_msd.py b/tests/test_pl_rdf_msd.py
new file mode 100644
index 0000000..cc2e2aa
--- /dev/null
+++ b/tests/test_pl_rdf_msd.py
@@ -0,0 +1,97 @@
+"""mv.pl.rdf_msd, and the MSD trace mv.md.run now keeps for it."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import matverse as mv
+
+
+@pytest.fixture(autouse=True)
+def _agg():
+    import matplotlib
+    matplotlib.use("Agg")
+
+
+def _frames(md, n=20):
+    cell = mv.structures(md, "input")[0]
+    frames = np.tile(np.array(cell.frac_coords), (n, 1, 1))
+    return frames + np.random.default_rng(0).normal(0, 0.01, frames.shape)
+
+
+@pytest.fixture(scope="module")
+def dynamic():
+    md = mv.datasets.metals(["Cu"], supercell=(2, 2, 2))
+    mv.pp.describe(md)
+    mv.md.run(md, level="emt", steps=40, equilibration=10, sample_every=5)
+    return md
+
+
+@pytest.fixture(scope="module")
+def diffusing(dynamic):
+    md = dynamic.copy()
+    pytest.importorskip("pymatgen.analysis.diffusion.aimd.rdf")
+    mv.md.rdf(md, _frames(md), species="Cu", reference="Cu", r_max=6.0)
+    return md
+
+
+class TestTheProducerKeepsTheTrace:
+    def test_run_deposits_the_msd_on_the_sampling_grid(self, dynamic):
+        assert "md_msd_trace_emt" in dynamic.obsm
+        trace = dynamic.obsm["md_msd_trace_emt"]
+        assert trace.shape == (1, 8)
+        assert np.isfinite(trace).all() and (trace >= 0).all()
+        np.testing.assert_array_equal(
+            mv.grid_of(dynamic, "md_msd_trace"),
+            mv.grid_of(dynamic, "md_temperature_trace"))
+        assert dynamic.uns["grids"]["md_msd_trace"]["unit"] == "ps"
+        assert trace[0, -1] == pytest.approx(
+            float(dynamic.obs["msd_emt"].iloc[0])), \
+            "the scalar MSD is the last point of the trace"
+
+    def test_a_rerun_of_a_different_length_replaces_the_trace(self, dynamic):
+        md = dynamic.copy()
+        mv.md.run(md, level="emt", steps=20, equilibration=10, sample_every=5)
+        assert md.obsm["md_msd_trace_emt"].shape == (1, 4)
+        assert len(mv.grid_of(md, "md_msd_trace")) == 4
+
+    def test_the_claim_is_in_the_registry(self):
+        entry = mv.registry.get("mv.md.run")
+        assert "md_msd_trace_{level}" in entry["produces"]["obsm"]
+
+
+class TestRdfMsd:
+    def test_both_panels_with_units(self, diffusing):
+        ax = mv.pl.rdf_msd(diffusing, level="emt")
+        msd = ax._matverse_msd_ax
+        assert ax.figure is msd.figure and ax is not msd
+        assert "Å" in ax.get_xlabel() and ax.get_ylabel() == "g(r)"
+        assert "ps" in msd.get_xlabel() and "Ų" in msd.get_ylabel()
+        assert len(ax.lines) >= 1
+        assert len(msd.lines) == 2, "the MSD and the 6Dt line"
+        assert "D = " in msd.get_legend().get_texts()[0].get_text()
+        assert "Cu" in ax.get_title() and "emt" in msd.get_title()
+
+    def test_the_slope_is_the_recorded_diffusivity(self, diffusing):
+        from matverse.md import _A2_PER_PS_TO_CM2_PER_S
+        msd = mv.pl.rdf_msd(diffusing, level="emt")._matverse_msd_ax
+        dashed = [l for l in msd.lines if l.get_linestyle() == "--"][0]
+        x, y = dashed.get_xdata(), dashed.get_ydata()
+        slope = (y[-1] - y[0]) / (x[-1] - x[0])
+        d = float(diffusing.obs["diffusivity_emt"].iloc[0])
+        assert slope == pytest.approx(6.0 * d / _A2_PER_PS_TO_CM2_PER_S)
+
+    def test_one_panel_alone(self, dynamic, diffusing):
+        ax = mv.pl.rdf_msd(dynamic, level="emt", which="msd")
+        assert "Ų" in ax.get_ylabel() and not hasattr(ax, "_matverse_msd_ax")
+        ax = mv.pl.rdf_msd(diffusing, which="rdf")
+        assert ax.get_ylabel() == "g(r)"
+
+    def test_says_which_producer_is_missing(self, dynamic, diffusing):
+        with pytest.raises(ValueError, match="mv.md.rdf"):
+            mv.pl.rdf_msd(dynamic, level="emt")
+        with pytest.raises(ValueError, match="mv.md.run"):
+            mv.pl.rdf_msd(diffusing, level="not-run")
+        with pytest.raises(ValueError, match="which must be"):
+            mv.pl.rdf_msd(diffusing, which="all")
diff --git a/tests/test_pl_wulff.py b/tests/test_pl_wulff.py
new file mode 100644
index 0000000..b080fd2
--- /dev/null
+++ b/tests/test_pl_wulff.py
@@ -0,0 +1,84 @@
+"""mv.pl.wulff, and the polyhedron mv.surf.wulff now keeps for it."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import matverse as mv
+
+
+@pytest.fixture(autouse=True)
+def _agg():
+    import matplotlib
+    matplotlib.use("Agg")
+
+
+@pytest.fixture(scope="module")
+def shaped():
+    bulk = mv.datasets.metals(["Cu"])
+    mv.pp.describe(bulk)
+    mv.calc.energy(bulk, level="emt")
+    facets = mv.surf.slabs(bulk, max_index=1)
+    mv.calc.energy(facets, level="emt")
+    mv.surf.surface_energy(facets, bulk, level="emt")
+    mv.surf.wulff(facets, bulk, level="emt")
+    return bulk, facets
+
+
+class TestTheProducerKeepsTheShape:
+    def test_every_face_on_the_particle_is_stored(self, shaped):
+        bulk, facets = shaped
+        shape = bulk.uns["wulff"]["emt"]["0"]
+        families = list(shape["face_miller"])
+        # fcc cut to max_index=1: 8 (111), 6 (100), 12 (110) faces.
+        assert sorted(set(families)) == ["1_0_0", "1_1_0", "1_1_1"]
+        assert len(families) == 26
+        vertices = np.asarray(shape["vertices"])
+        index = np.asarray(shape["face_index"])
+        assert vertices.shape == (len(index), 3)
+        assert index.max() == 25 and (np.bincount(index) >= 3).all()
+        expressed = facets.obs.set_index("miller")["wulff_area_fraction_emt"]
+        for miller, fraction in shape["area_fractions"].items():
+            assert fraction == pytest.approx(float(expressed[miller]))
+
+    def test_the_claim_is_back_in_the_registry(self):
+        entry = mv.registry.get("mv.surf.wulff")
+        assert entry["produces"]["bulk.uns"] == ["wulff"]
+
+    def test_it_survives_h5ad(self, shaped, tmp_path):
+        import anndata
+        bulk, _ = shaped
+        bulk.write_h5ad(tmp_path / "bulk.h5ad")
+        back = anndata.read_h5ad(tmp_path / "bulk.h5ad")
+        assert np.asarray(back.uns["wulff"]["emt"]["0"]["vertices"]).shape \
+            == np.asarray(bulk.uns["wulff"]["emt"]["0"]["vertices"]).shape
+
+
+class TestWulff:
+    def test_draws_the_polyhedron_coloured_by_family(self, shaped):
+        bulk, _ = shaped
+        ax = mv.pl.wulff(bulk, level="emt")
+        assert ax.name == "3d"
+        assert ax._matverse_n_faces == 26
+        (collection,) = ax.collections
+        assert len(collection.get_facecolor()) == 26
+        labels = [t.get_text() for t in ax.get_legend().get_texts()]
+        assert len(labels) == 3
+        assert any(l.startswith("(1 1 1)") and "% of the surface" in l
+                   for l in labels)
+        assert "Cu" in ax.get_title() and "emt" in ax.get_title()
+
+    def test_a_row_can_be_named(self, shaped):
+        bulk, _ = shaped
+        assert mv.pl.wulff(bulk, level="emt", row="Cu")._matverse_n_faces == 26
+
+    def test_says_what_is_missing(self, shaped):
+        bulk, _ = shaped
+        bare = mv.datasets.metals(["Cu"])
+        with pytest.raises(ValueError, match="mv.surf.wulff"):
+            mv.pl.wulff(bare, level="emt")
+        stale = bulk.copy()
+        stale.uns["wulff"]["emt"]["0"] = {"n_facets": 3, "anisotropy": 0.05}
+        with pytest.raises(ValueError, match="older mv.surf.wulff"):
+            mv.pl.wulff(stale, level="emt")