diff --git a/src/py4vasp/_calculation/__init__.py b/src/py4vasp/_calculation/__init__.py
index 0fdb11e8..dc6d3e33 100644
--- a/src/py4vasp/_calculation/__init__.py
+++ b/src/py4vasp/_calculation/__init__.py
@@ -380,12 +380,10 @@ def __dir__(self):
     def _compute_database_data(self) -> dict:
         """Iterate over all quantities in _REGISTRY and collect database properties.
 
-        Returns a flat dict mapping property keys to handler result dicts.
-        Keys use the format:
-        - ``<quantity>`` for the default selection
-        - ``<quantity>_<selection>`` for non-default selections
-        - ``<group>_<quantity>`` / ``<group>_<quantity>_<selection>`` for groups
-        Leading underscores are stripped from private quantity names.
+        Returns a nested dict ``{quantity: {selection: model}}``. The outer key is
+        the (underscore-stripped) quantity name; the inner dict is keyed by
+        selection with the default source keyed ``"default"``. Group members use
+        ``<group>_<quantity>`` (e.g. ``phonon_mode``) as their outer key.
         """
         _ensure_all_quantities_imported()
         properties = {}
@@ -435,11 +433,12 @@ def _ensure_all_quantities_imported():
 
 
 def _collect_to_database(dispatcher_cls, source, properties):
-    """Call dispatcher._to_database() and merge results into *properties*.
+    """Call dispatcher._to_database() and deep-merge results into *properties*.
 
-    Each dispatcher already keys its results by the full quantity name; group members
-    use ``<group>_<member>`` via their ``_quantity_name`` (e.g. ``phonon_mode``), so the
-    results merge directly without any additional group prefix.
+    Each dispatcher returns a nested ``{quantity: {selection: model}}`` dict keyed by
+    the full quantity name; group members use ``<group>_<member>`` via their
+    ``_quantity_name`` (e.g. ``phonon_mode``). The nested dicts are merged so that
+    quantities and their selections accumulate without clobbering each other.
     """
     dispatcher = dispatcher_cls(
         source=source, quantity_name=dispatcher_cls._quantity_name
@@ -452,7 +451,8 @@ def _collect_to_database(dispatcher_cls, source, properties):
         return
     except Exception:
         return
-    properties.update(result)
+    for quantity, selections in result.items():
+        properties.setdefault(quantity, {}).update(selections)
 
 
 def _rebuild_public_registry_views():
diff --git a/src/py4vasp/_calculation/_dispersion.py b/src/py4vasp/_calculation/_dispersion.py
index 3c792814..9e86aef0 100644
--- a/src/py4vasp/_calculation/_dispersion.py
+++ b/src/py4vasp/_calculation/_dispersion.py
@@ -8,7 +8,6 @@
     DataSource,
     merge_default,
     merge_strings,
-    merge_to_database,
     quantity,
 )
 from py4vasp._calculation.kpoint import KpointHandler
@@ -146,14 +145,9 @@ def plot(self, projections=None):
             projections,
         )
 
-    def _to_database(self) -> dict:
-        """Return {quantity[_selection]: handler_result} for database storage."""
-        return merge_to_database(
-            self._source,
-            self._quantity_name,
-            DispersionHandler.from_data,
-            DispersionHandler.to_database,
-        )
+    # Dispersion has no standalone database entry: its data is folded into the band
+    # and phonon_band models (see BandHandler / PhononBandHandler.to_database), so the
+    # dispatcher deliberately exposes no `_to_database`.
 
 
 def _band_structure(data, projections):
diff --git a/src/py4vasp/_calculation/_stoichiometry.py b/src/py4vasp/_calculation/_stoichiometry.py
index 726ca164..fc0d8fa8 100644
--- a/src/py4vasp/_calculation/_stoichiometry.py
+++ b/src/py4vasp/_calculation/_stoichiometry.py
@@ -9,7 +9,6 @@
     DataSource,
     merge_default,
     merge_strings,
-    merge_to_database,
     quantity,
 )
 from py4vasp._calculation.selection import Selection
@@ -417,14 +416,9 @@ def number_atoms(self):
             StoichiometryHandler.number_atoms,
         )
 
-    def _to_database(self) -> dict:
-        """Return {quantity[_selection]: handler_result} for database storage."""
-        return merge_to_database(
-            self._source,
-            self._quantity_name,
-            StoichiometryHandler.from_data,
-            StoichiometryHandler.to_database,
-        )
+    # Stoichiometry has no standalone database entry: its data is folded into the
+    # structure model (see StructureHandler.to_database), so the dispatcher
+    # deliberately exposes no `_to_database`.
 
 
 def raw_stoichiometry_from_ase(structure):
diff --git a/src/py4vasp/_calculation/band.py b/src/py4vasp/_calculation/band.py
index 85afbab4..eecab7cf 100644
--- a/src/py4vasp/_calculation/band.py
+++ b/src/py4vasp/_calculation/band.py
@@ -106,6 +106,8 @@ def to_database(self, selection=None, fermi_energy=None) -> Band_DB:
             else None
         )
 
+        dispersion = self._dispersion().to_database()
+
         return Band_DB(
             num_considered_bands=num_checked_bands,
             num_occupied_bands=num_total_occupied,
@@ -113,6 +115,12 @@ def to_database(self, selection=None, fermi_energy=None) -> Band_DB:
             num_occupied_bands_down=num_occupied_down,
             fermi_energy_raw=raw_fermi_energy,
             fermi_energy=fermi_energy or raw_fermi_energy,
+            eigenvalue_min=dispersion.eigenvalue_min,
+            eigenvalue_max=dispersion.eigenvalue_max,
+            eigenvalue_min_up=dispersion.eigenvalue_min_up,
+            eigenvalue_max_up=dispersion.eigenvalue_max_up,
+            eigenvalue_min_down=dispersion.eigenvalue_min_down,
+            eigenvalue_max_down=dispersion.eigenvalue_max_down,
         )
 
     def to_graph(self, selection=None, fermi_energy=None, width=0.5) -> graph.Graph:
diff --git a/src/py4vasp/_calculation/current_density.py b/src/py4vasp/_calculation/current_density.py
index cc009cc0..d2df5303 100644
--- a/src/py4vasp/_calculation/current_density.py
+++ b/src/py4vasp/_calculation/current_density.py
@@ -16,8 +16,9 @@
 )
 from py4vasp._calculation.structure import StructureHandler
 from py4vasp._raw import data as raw
+from py4vasp._raw.data_db import CurrentDensity_DB
 from py4vasp._third_party import graph
-from py4vasp._util import documentation, import_, slicing
+from py4vasp._util import check, documentation, import_, slicing
 from py4vasp._util.density import SliceArguments, Visualizer
 
 pretty = import_.optional("IPython.lib.pretty")
@@ -62,8 +63,29 @@ def to_dict(self) -> dict:
             **self._read_current_densities(),
         }
 
-    def to_database(self) -> dict:
-        return {}
+    def to_database(self) -> "CurrentDensity_DB":
+        """Serialize a scalar summary of the current density for database storage."""
+        grid_shape, magnitudes = None, []
+        for key in self._raw_current_density.valid_indices:
+            current_density = self._raw_current_density[key].current_density
+            if check.is_none(current_density):
+                continue
+            # current_density has shape (3, nz, ny, nx); |j| is the norm over the
+            # vector component axis, leaving the grid.
+            vector_field = np.array(current_density[:])
+            magnitudes.append(np.linalg.norm(vector_field, axis=0).ravel())
+            if grid_shape is None:
+                nz, ny, nx = vector_field.shape[1:]
+                grid_shape = [nx, ny, nz]
+        if not magnitudes:
+            return CurrentDensity_DB()
+        magnitude = np.concatenate(magnitudes)
+        return CurrentDensity_DB(
+            grid_shape=grid_shape,
+            magnitude_min=float(np.min(magnitude)),
+            magnitude_max=float(np.max(magnitude)),
+            magnitude_mean=float(np.mean(magnitude)),
+        )
 
     def to_contour(
         self,
diff --git a/src/py4vasp/_calculation/dispatch.py b/src/py4vasp/_calculation/dispatch.py
index 82cc527e..277ca948 100644
--- a/src/py4vasp/_calculation/dispatch.py
+++ b/src/py4vasp/_calculation/dispatch.py
@@ -389,9 +389,10 @@ def merge_to_database(
 
     Returns
     -------
-    dict[str, result]
-        Maps ``quantity_name`` (default source) or ``quantity_name_source``
-        (other sources) to each handler result.
+    dict[str, dict[str, result]]
+        Nested ``{quantity: {selection: result}}``. The outer key is the
+        (underscore-stripped) quantity name; the inner dict is keyed by selection
+        with the default source keyed ``"default"``. Empty quantities are omitted.
     """
     wrapped_source = SuppressErrorsSourceWrapper(source)
     base = (key_name or quantity_name).lstrip("_")
@@ -405,12 +406,11 @@ def merge_to_database(
         *args,
         **kwargs,
     )
-    results = {
-        base if sel == "default" else f"{base}_{sel}": result
-        for sel, result in raw_results.items()
-        if _result_has_data(result)
+    selections = {
+        sel: result for sel, result in raw_results.items() if _result_has_data(result)
     }
-    return _drop_duplicates_of_default(results, base)
+    selections = _drop_duplicates_of_default(selections)
+    return {base: selections} if selections else {}
 
 
 def _all_sources_selection(quantity_name):
@@ -426,15 +426,15 @@ def _all_sources_selection(quantity_name):
     return ", ".join(str(source_name) for source_name in sources)
 
 
-def _drop_duplicates_of_default(results, base):
-    """Drop non-default entries whose result equals the default result."""
-    if base not in results:
-        return results
-    default_result = results[base]
+def _drop_duplicates_of_default(selections):
+    """Drop non-default selections whose result equals the ``default`` result."""
+    if "default" not in selections:
+        return selections
+    default_result = selections["default"]
     return {
-        key: result
-        for key, result in results.items()
-        if key == base or result != default_result
+        sel: result
+        for sel, result in selections.items()
+        if sel == "default" or result != default_result
     }
 
 
diff --git a/src/py4vasp/_calculation/pair_correlation.py b/src/py4vasp/_calculation/pair_correlation.py
index ffdf2911..da7e622a 100644
--- a/src/py4vasp/_calculation/pair_correlation.py
+++ b/src/py4vasp/_calculation/pair_correlation.py
@@ -2,6 +2,8 @@
 # Licensed under the Apache License 2.0 (http://www.apache.org/licenses/LICENSE-2.0)
 import copy
 
+import numpy as np
+
 from py4vasp import raw
 from py4vasp._calculation import slice_
 from py4vasp._calculation.dispatch import (
@@ -17,6 +19,11 @@
 from py4vasp._third_party import graph
 from py4vasp._util import check, convert, documentation, index, select
 
+# Minimum height of the total pair correlation function for a local maximum to be
+# considered the first peak. The value corresponds to the ideal-gas baseline, so
+# small bumps below it (e.g. numerical noise before the first shell) are ignored.
+_FIRST_PEAK_THRESHOLD = 1.0
+
 
 def _selection_string(default):
     return f"""\
@@ -62,13 +69,42 @@ def labels(self) -> tuple:
         """Return all possible labels for the selection string."""
         return tuple(convert.text_to_string(label) for label in self._raw_data.labels)
 
-    def to_database(self) -> dict:
+    def to_database(self) -> PairCorrelation_DB:
         """Serialize pair-correlation data for database storage."""
         distance_min, distance_max = None, None
         if not check.is_none(self._raw_data.distances):
             distance_min = float(self._raw_data.distances[0])
             distance_max = float(self._raw_data.distances[-1])
-        return PairCorrelation_DB(distance_min=distance_min, distance_max=distance_max)
+        first_peak_position, first_peak_height = self._first_peak()
+        return PairCorrelation_DB(
+            distance_min=distance_min,
+            distance_max=distance_max,
+            first_peak_position=first_peak_position,
+            first_peak_height=first_peak_height,
+        )
+
+    def _first_peak(self):
+        """Position and height of the first peak of the total pair correlation.
+
+        Scans the existing sample points for the first strict local maximum of the
+        total g(r) whose height exceeds :data:`_FIRST_PEAK_THRESHOLD`, so small
+        bumps before the first real shell are skipped. Returns ``(None, None)`` when
+        no such peak exists.
+        """
+        if check.is_none(self._raw_data.distances) or check.is_none(
+            self._raw_data.function
+        ):
+            return None, None
+        distances = np.asarray(self._raw_data.distances[:])
+        # index 0 of the function is the total pair correlation by convention
+        total = np.asarray(self._raw_data.function)[self._steps_or_last, 0]
+        if distances.shape != total.shape or total.size < 3:
+            return None, None
+        for i in range(1, total.size - 1):
+            is_local_max = total[i] > total[i - 1] and total[i] > total[i + 1]
+            if is_local_max and total[i] > _FIRST_PEAK_THRESHOLD:
+                return float(distances[i]), float(total[i])
+        return None, None
 
     @property
     def _steps_or_last(self):
diff --git a/src/py4vasp/_calculation/phonon_band.py b/src/py4vasp/_calculation/phonon_band.py
index c7888911..ca222a80 100644
--- a/src/py4vasp/_calculation/phonon_band.py
+++ b/src/py4vasp/_calculation/phonon_band.py
@@ -17,6 +17,7 @@
     merge_to_database,
     quantity,
 )
+from py4vasp._raw.data_db import PhononBand_DB
 from py4vasp._third_party import graph
 from py4vasp._util import convert, documentation, index, select
 
@@ -46,8 +47,12 @@ def to_dict(self) -> dict:
             "modes": self._modes(),
         }
 
-    def to_database(self) -> dict:
-        return {}
+    def to_database(self) -> PhononBand_DB:
+        dispersion = self._dispersion().to_database()
+        return PhononBand_DB(
+            eigenvalue_min=dispersion.eigenvalue_min,
+            eigenvalue_max=dispersion.eigenvalue_max,
+        )
 
     def to_graph(self, selection=None, width=1.0) -> graph.Graph:
         projections = self._projections(selection, width)
diff --git a/src/py4vasp/_calculation/structure.py b/src/py4vasp/_calculation/structure.py
index f3424478..8990ce3c 100644
--- a/src/py4vasp/_calculation/structure.py
+++ b/src/py4vasp/_calculation/structure.py
@@ -13,12 +13,14 @@
 from py4vasp._calculation.cell import CellHandler
 from py4vasp._calculation.dispatch import (
     DataSource,
+    SuppressErrorsSourceWrapper,
+    _result_has_data,
     merge_default,
     merge_strings,
-    merge_to_database,
     quantity,
 )
-from py4vasp._raw.data_db import Structure_DB
+from py4vasp._raw.data_db import Stoichiometry_DB, Structure_DB
+from py4vasp._raw.definition import unique_selections as _schema_unique_selections
 from py4vasp._third_party import view
 from py4vasp._util import check, import_, parse
 
@@ -35,6 +37,9 @@
     TypeError,
     ValueError,
     IndexError,
+    # a source may reference datasets that are absent/malformed in the HDF5 file;
+    # reading them can raise a low-level RuntimeError which we treat as "no data"
+    RuntimeError,
 )
 
 
@@ -220,148 +225,94 @@ def number_steps(self) -> int:
         else:
             return 1
 
-    def to_database(self) -> dict:
-        """Return database-ready structure data."""
-        # Temporarily use all steps for database
+    def to_database(self, steps=-1) -> Structure_DB:
+        """Return database-ready data for a single structure geometry.
+
+        *steps* selects which geometry to describe (default ``-1``, the final step).
+        The database splits a calculation into separate ``initial`` and ``final``
+        structure models, so each :class:`Structure_DB` holds one geometry with
+        unprefixed fields.
+        """
+        # Temporarily override the steps used for the database
         saved_steps = self._steps
         saved_is_slice = self._is_slice
         saved_slice = self._slice
-        self._steps = slice(None)
-        self._is_slice = True
-        self._slice = slice(None)
+        self._steps = steps
+        self._is_slice = isinstance(steps, slice)
+        if self._is_slice:
+            self._slice = steps
+        elif steps == -1:
+            self._slice = slice(-1, None)
+        else:
+            self._slice = slice(steps, steps + 1)
 
-        final_lattice, initial_lattice = ([None, None, None] for _ in range(2))
+        try:
+            return self._single_geometry_database()
+        finally:
+            self._steps = saved_steps
+            self._is_slice = saved_is_slice
+            self._slice = saved_slice
+
+    def _single_geometry_database(self) -> Structure_DB:
+        lattice = [None, None, None]
         with suppress(*_TO_DATABASE_SUPPRESSED_EXCEPTIONS):
             lattices = self.lattice_vectors()
-            final_lattice = lattices[-1] if lattices.ndim == 3 else lattices
-            initial_lattice = lattices[0] if lattices.ndim == 3 else lattices
-            if final_lattice.ndim != 2:
-                final_lattice = [None, None, None]
-            if initial_lattice.ndim != 2:
-                initial_lattice = [None, None, None]
-
-        volume_final, volume_initial = (None for _ in range(2))
+            lattice = lattices[-1] if lattices.ndim == 3 else lattices
+            if lattice.ndim != 2:
+                lattice = [None, None, None]
+
+        volume = None
         with suppress(*_TO_DATABASE_SUPPRESSED_EXCEPTIONS):
             volumes = self.volume()
-            volume_final = (
+            volume = (
                 volumes[-1]
                 if not isinstance(volumes, (float, np.float64, np.float32))
                 else volumes
             )
-            volume_initial = (
-                volumes[0]
-                if not isinstance(volumes, (float, np.float64, np.float32))
-                else volumes
-            )
 
-        lengths_final, angles_final, lengths_initial, angles_initial = (
-            None for _ in range(4)
-        )
-        (
-            cell_area_2d_final,
-            cell_area_2d_span_final,
-            cell_area_2d_initial,
-            cell_area_2d_span_initial,
-        ) = (None for _ in range(4))
+        lengths = angles = None
+        cell_area_2d = cell_area_2d_span = None
         dimensionality = 3
         with suppress(*_TO_DATABASE_SUPPRESSED_EXCEPTIONS):
             dimensionality = self._dimensionality()
 
         with suppress(*_TO_DATABASE_SUPPRESSED_EXCEPTIONS):
             cell_ = self._cell()
-            lengths = cell_.lengths()
-            lengths_final = lengths[-1] if lengths.ndim == 2 else lengths
-            lengths_initial = lengths[0] if lengths.ndim == 2 else lengths
-            angles = cell_.angles()
-            angles_final = angles[-1] if angles.ndim == 2 else angles
-            angles_initial = angles[0] if angles.ndim == 2 else angles
+            all_lengths = cell_.lengths()
+            lengths = all_lengths[-1] if all_lengths.ndim == 2 else all_lengths
+            all_angles = cell_.angles()
+            angles = all_angles[-1] if all_angles.ndim == 2 else all_angles
             if dimensionality == 2:
-                cell_area_2d, cell_area_2d_span = cell_._area_2d()
-                cell_area_2d_final = (
-                    cell_area_2d[-1]
-                    if isinstance(cell_area_2d, np.ndarray)
-                    else cell_area_2d
-                )
-                cell_area_2d_initial = (
-                    cell_area_2d[0]
-                    if isinstance(cell_area_2d, np.ndarray)
-                    else cell_area_2d
-                )
-                cell_area_2d_span_final = (
-                    cell_area_2d_span[-1]
-                    if isinstance(cell_area_2d_span, list)
-                    else cell_area_2d_span
-                )
-                cell_area_2d_span_initial = (
-                    cell_area_2d_span[0]
-                    if isinstance(cell_area_2d_span, list)
-                    else cell_area_2d_span
-                )
+                area, span = cell_._area_2d()
+                cell_area_2d = area[-1] if isinstance(area, np.ndarray) else area
+                cell_area_2d_span = span[-1] if isinstance(span, list) else span
 
         num_atoms = self.number_atoms() or None
 
-        # Restore steps
-        self._steps = saved_steps
-        self._is_slice = saved_is_slice
-        self._slice = saved_slice
+        stoichiometry = Stoichiometry_DB()
+        with suppress(*_TO_DATABASE_SUPPRESSED_EXCEPTIONS):
+            stoichiometry = self._stoichiometry().to_database()
 
         return Structure_DB(
             num_ions=num_atoms,
             dimensionality=dimensionality,
-            final_cell_volume=volume_final,
-            final_cell_area_2d=cell_area_2d_final,
-            final_cell_area_2d_span=cell_area_2d_span_final,
-            final_lattice_vector_1=(
-                list(final_lattice[0]) if final_lattice[0] is not None else None
-            ),
-            final_lattice_vector_2=(
-                list(final_lattice[1]) if final_lattice[1] is not None else None
-            ),
-            final_lattice_vector_3=(
-                list(final_lattice[2]) if final_lattice[2] is not None else None
-            ),
-            final_lattice_vector_1_length=(
-                lengths_final[0] if lengths_final is not None else None
-            ),
-            final_lattice_vector_2_length=(
-                lengths_final[1] if lengths_final is not None else None
-            ),
-            final_lattice_vector_3_length=(
-                lengths_final[2] if lengths_final is not None else None
-            ),
-            final_angle_alpha=(angles_final[0] if angles_final is not None else None),
-            final_angle_beta=(angles_final[1] if angles_final is not None else None),
-            final_angle_gamma=(angles_final[2] if angles_final is not None else None),
-            initial_cell_volume=volume_initial,
-            initial_cell_area_2d=cell_area_2d_initial,
-            initial_cell_area_2d_span=cell_area_2d_span_initial,
-            initial_lattice_vector_1=(
-                list(initial_lattice[0]) if initial_lattice[0] is not None else None
-            ),
-            initial_lattice_vector_2=(
-                list(initial_lattice[1]) if initial_lattice[1] is not None else None
-            ),
-            initial_lattice_vector_3=(
-                list(initial_lattice[2]) if initial_lattice[2] is not None else None
-            ),
-            initial_lattice_vector_1_length=(
-                lengths_initial[0] if lengths_initial is not None else None
-            ),
-            initial_lattice_vector_2_length=(
-                lengths_initial[1] if lengths_initial is not None else None
-            ),
-            initial_lattice_vector_3_length=(
-                lengths_initial[2] if lengths_initial is not None else None
-            ),
-            initial_angle_alpha=(
-                angles_initial[0] if angles_initial is not None else None
-            ),
-            initial_angle_beta=(
-                angles_initial[1] if angles_initial is not None else None
-            ),
-            initial_angle_gamma=(
-                angles_initial[2] if angles_initial is not None else None
-            ),
+            ion_types=stoichiometry.ion_types,
+            num_ion_types=stoichiometry.num_ion_types,
+            num_ion_types_primitive=stoichiometry.num_ion_types_primitive,
+            formula=stoichiometry.formula,
+            compound=stoichiometry.compound,
+            cell_volume=volume,
+            cell_area_2d=cell_area_2d,
+            cell_area_2d_span=cell_area_2d_span,
+            lattice_vector_1=list(lattice[0]) if lattice[0] is not None else None,
+            lattice_vector_2=list(lattice[1]) if lattice[1] is not None else None,
+            lattice_vector_3=list(lattice[2]) if lattice[2] is not None else None,
+            lattice_vector_1_length=lengths[0] if lengths is not None else None,
+            lattice_vector_2_length=lengths[1] if lengths is not None else None,
+            lattice_vector_3_length=lengths[2] if lengths is not None else None,
+            angle_alpha=angles[0] if angles is not None else None,
+            angle_beta=angles[1] if angles is not None else None,
+            angle_gamma=angles[2] if angles is not None else None,
         )
 
     def _dimensionality(self) -> Union[int, np.ndarray]:
@@ -1229,13 +1180,98 @@ def number_steps(self):
         )
 
     def _to_database(self) -> dict:
-        """Return {quantity[_selection]: handler_result} for database storage."""
-        return merge_to_database(
-            self._source,
-            self._quantity_name,
-            StructureHandler.from_data,
-            StructureHandler.to_database,
-        )
+        """Return ``{"structure": {selection: Structure_DB}}`` for the database.
+
+        Uses a custom merge instead of the generic :func:`merge_to_database`: the
+        final structure is always stored, the initial structure (first ionic step
+        of the ``default`` trajectory) only when it differs from the final one, and
+        any additional source only when it differs from both.
+        """
+        models = _collect_structure_database(self._source, self._quantity_name)
+        return {"structure": models} if models else {}
+
+
+def _collect_structure_database(source, quantity_name):
+    """Build ``{selection: Structure_DB}`` from the available structure sources.
+
+    The ``final`` source (or, when absent, the last step of the ``default``
+    trajectory) provides the final structure, which is always included. The first
+    ionic step of the ``default`` trajectory provides the initial structure, kept
+    only when it differs from the final one. Every other source is kept only when
+    it differs from both.
+    """
+    entries = _read_structure_entries(source, quantity_name)
+    return _assemble_structure_models(entries)
+
+
+def _read_structure_entries(source, quantity_name):
+    """Materialize ``{source: {step: (model, geometry)}}`` for every schema source.
+
+    The reads happen while the source's access context is open (raw data holds
+    lazy references that become invalid once the file is closed) and the arrays
+    are copied so they survive after the context exits. The ``default`` trajectory
+    is read at both its first (``0``) and last (``-1``) step; every other source
+    only at its final step (``-1``).
+    """
+    wrapped = SuppressErrorsSourceWrapper(source)
+    entries = {}
+    for name in _schema_unique_selections(quantity_name.lstrip("_")):
+        selection = None if name == "default" else name
+        steps = (-1, 0) if name == "default" else (-1,)
+        with wrapped.access(quantity_name, selection=selection) as raw_structure:
+            if raw_structure is None:
+                continue
+            per_step = {}
+            for step in steps:
+                entry = _build_structure_entry(raw_structure, step)
+                if entry is not None:
+                    per_step[step] = entry
+            if per_step:
+                entries[name] = per_step
+    return entries
+
+
+def _build_structure_entry(raw_structure, step):
+    """Return ``(Structure_DB, (lattice, positions))`` for a single step, or None."""
+    with suppress(*_TO_DATABASE_SUPPRESSED_EXCEPTIONS):
+        handler = StructureHandler.from_data(raw_structure, steps=step)
+        geometry = (np.array(handler.lattice_vectors()), np.array(handler.positions()))
+        model = StructureHandler.from_data(raw_structure).to_database(steps=step)
+        if not _result_has_data(model):
+            return None
+        return model, geometry
+    return None
+
+
+def _assemble_structure_models(entries):
+    """Apply the final/initial/other rules to the collected per-source entries."""
+    models, geometries = {}, {}
+    final = entries.get("final", {}).get(-1) or entries.get("default", {}).get(-1)
+    if final is not None:
+        models["final"], geometries["final"] = final
+    initial = entries.get("default", {}).get(0)
+    if initial is not None and not _duplicate_geometry(initial[1], geometries):
+        models["initial"], geometries["initial"] = initial
+    for name, per_step in entries.items():
+        if name in ("default", "final") or -1 not in per_step:
+            continue
+        model, geometry = per_step[-1]
+        if _duplicate_geometry(geometry, geometries):
+            continue
+        models[name], geometries[name] = model, geometry
+    return models
+
+
+def _duplicate_geometry(geometry, geometries):
+    """Whether *geometry* matches any already-collected geometry."""
+    return any(_same_geometry(geometry, other) for other in geometries.values())
+
+
+def _same_geometry(first, second):
+    """Two geometries are equal when both lattice vectors and positions match."""
+    return all(
+        np.shape(a) == np.shape(b) and np.allclose(a, b) for a, b in zip(first, second)
+    )
 
 
 def _cell_from_ase(structure):
diff --git a/src/py4vasp/_raw/data.py b/src/py4vasp/_raw/data.py
index 86ee439e..9e1c2dd3 100644
--- a/src/py4vasp/_raw/data.py
+++ b/src/py4vasp/_raw/data.py
@@ -185,10 +185,10 @@ class _DatabaseData:
     metadata: CalculationMetaData
     properties: dict = dataclasses.field(default_factory=dict)
     """Properties extracted from the calculation.
-    Keys follow the format 'quantity' (default selection) or 'quantity_selection'
-    (non-default selection). Group quantities use 'group_quantity' or
-    'group_quantity_selection'. Leading underscores are stripped from private
-    quantity names."""
+    A nested dict ``{quantity: {selection: model}}``. The outer key is the quantity
+    name (group members use 'group_quantity', e.g. 'phonon_mode'); leading
+    underscores are stripped from private quantity names. The inner dict is keyed by
+    selection with the default source keyed 'default'."""
 
 
 @dataclasses.dataclass
diff --git a/src/py4vasp/_raw/data_db.py b/src/py4vasp/_raw/data_db.py
index ce27081b..1ccfe088 100644
--- a/src/py4vasp/_raw/data_db.py
+++ b/src/py4vasp/_raw/data_db.py
@@ -78,6 +78,20 @@ class Band_DB(_DBDataMixin):
     fermi_energy: Optional[float] = None
     """The Fermi energy used for plotting the band structure, which may be different from the raw Fermi energy based on user input."""
 
+    # dispersion data folded into the band model
+    eigenvalue_min: Optional[float] = None
+    """The minimum eigenvalue across all bands and k-points, in eV."""
+    eigenvalue_max: Optional[float] = None
+    """The maximum eigenvalue across all bands and k-points, in eV."""
+    eigenvalue_min_up: Optional[float] = None
+    """The minimum eigenvalue for spin-up electrons across all bands and k-points, in eV."""
+    eigenvalue_max_up: Optional[float] = None
+    """The maximum eigenvalue for spin-up electrons across all bands and k-points, in eV."""
+    eigenvalue_min_down: Optional[float] = None
+    """The minimum eigenvalue for spin-down electrons across all bands and k-points, in eV."""
+    eigenvalue_max_down: Optional[float] = None
+    """The maximum eigenvalue for spin-down electrons across all bands and k-points, in eV."""
+
 
 @dataclass
 class Bandgap_DB(_DBDataMixin):
@@ -164,6 +178,23 @@ class BornEffectiveCharge_DB(_DBDataMixin):
     """The index of the ion with the maximum eigenvalue of the Born effective charge tensor."""
 
 
+@dataclass
+class CurrentDensity_DB(_DBDataMixin):
+    """Data class for storing current density data in the database.
+
+    The current density is a vector field on a grid; the summaries below refer to its
+    magnitude |j| aggregated over all grid points and perturbations."""
+
+    grid_shape: Optional[List[int]] = field(default_factory=lambda: None)
+    """The shape of the grid on which the current density is evaluated, in the order (nx, ny, nz)."""
+    magnitude_min: Optional[float] = None
+    """The minimum magnitude of the current density across all grid points and perturbations."""
+    magnitude_max: Optional[float] = None
+    """The maximum magnitude of the current density across all grid points and perturbations."""
+    magnitude_mean: Optional[float] = None
+    """The mean magnitude of the current density across all grid points and perturbations."""
+
+
 @dataclass
 class DielectricFunction_DB(_DBDataMixin):
     """Data class for storing dielectric function data in the database."""
@@ -588,6 +619,10 @@ class PairCorrelation_DB(_DBDataMixin):
     """The minimum distance at which the pair correlation function was evaluated, in Å."""
     distance_max: Optional[float] = None
     """The maximum distance at which the pair correlation function was evaluated, in Å."""
+    first_peak_position: Optional[float] = None
+    """The distance of the first peak of the total pair correlation function, in Å."""
+    first_peak_height: Optional[float] = None
+    """The height (value of the total pair correlation function) at the first peak."""
 
 
 @dataclass
@@ -600,6 +635,18 @@ class PhononDos_DB(_DBDataMixin):
     """The maximum energy at which the phonon density of states was evaluated, in THz."""
 
 
+@dataclass
+class PhononBand_DB(_DBDataMixin):
+    """Data class for storing phonon band structure data in the database.
+
+    The dispersion (phonon frequencies) is folded into this model."""
+
+    eigenvalue_min: Optional[float] = None
+    """The minimum phonon frequency across all modes and q-points, in THz."""
+    eigenvalue_max: Optional[float] = None
+    """The maximum phonon frequency across all modes and q-points, in THz."""
+
+
 @dataclass
 class PhononMode_DB(_DBDataMixin):
     """Data class for storing phonon mode data in the database."""
@@ -830,73 +877,59 @@ class Stress_DB(_DBDataMixin):
 
 @dataclass
 class Structure_DB(_DBDataMixin):
-    """Data class for storing structure data in the database."""
+    """Data class for storing a single structure geometry in the database.
+
+    Each instance describes one geometry. The database entry for a calculation stores
+    a ``final`` structure and, when it differs, a separate ``initial`` structure rather
+    than packing both into one model, so the fields carry no ``initial_``/``final_``
+    prefix."""
 
     num_ions: Optional[int] = None
     """The number of ions in the structure."""
     dimensionality: Optional[int] = None
     """The dimensionality of the structure, as determined by the presence of vacuum along different lattice vectors.
-    
+
     - 3 = bulk structure with no vacuum along any lattice vector
     - 2 = slab structure with vacuum along one lattice vector
     - 1 = multi-atom molecule or wire structure with vacuum along two lattice vectors
     - 0 = single-atom structure with vacuum along all three lattice vectors"""
 
-    final_cell_volume: Optional[float] = None
-    """The volume of the unit cell on the final step, in Å^3."""
-    final_cell_area_2d: Optional[float] = None
+    cell_volume: Optional[float] = None
+    """The volume of the unit cell, in Å^3."""
+    cell_area_2d: Optional[float] = None
     """The area of the unit cell in 2D materials, calculated as the product of the two lattice vectors that are not along the vacuum direction, in Å^2."""
-    final_cell_area_2d_span: Optional[str] = None
+    cell_area_2d_span: Optional[str] = None
     """The two lattice vectors that are used to calculate the area of the unit cell in 2D materials, in the format '12', '13', or '23'."""
-    final_lattice_vector_1: Optional[List[float]] = field(default_factory=lambda: None)
-    """The first lattice vector on the final step, in Å."""
-    final_lattice_vector_2: Optional[List[float]] = field(default_factory=lambda: None)
-    """The second lattice vector on the final step, in Å."""
-    final_lattice_vector_3: Optional[List[float]] = field(default_factory=lambda: None)
-    """The third lattice vector on the final step, in Å."""
-    final_lattice_vector_1_length: Optional[float] = None
-    """The length of the first lattice vector on the final step, in Å."""
-    final_lattice_vector_2_length: Optional[float] = None
-    """The length of the second lattice vector on the final step, in Å."""
-    final_lattice_vector_3_length: Optional[float] = None
-    """The length of the third lattice vector on the final step, in Å."""
-    final_angle_alpha: Optional[float] = None
-    """The angle between the second and third lattice vectors on the final step, in degrees."""
-    final_angle_beta: Optional[float] = None
-    """The angle between the first and third lattice vectors on the final step, in degrees."""
-    final_angle_gamma: Optional[float] = None
-    """The angle between the first and second lattice vectors on the final step, in degrees."""
-
-    initial_cell_volume: Optional[float] = None
-    """The volume of the unit cell on the initial step, in Å^3."""
-    initial_cell_area_2d: Optional[float] = None
-    """The area of the unit cell in 2D materials on the initial step, calculated as the product of the two lattice vectors that are not along the vacuum direction, in Å^2."""
-    initial_cell_area_2d_span: Optional[str] = None
-    """The two lattice vectors that are used to calculate the area of the unit cell in 2D materials on the initial step, in the format '12', '13', or '23'."""
-    initial_lattice_vector_1: Optional[List[float]] = field(
-        default_factory=lambda: None
-    )
-    """The first lattice vector on the initial step, in Å."""
-    initial_lattice_vector_2: Optional[List[float]] = field(
-        default_factory=lambda: None
-    )
-    """The second lattice vector on the initial step, in Å."""
-    initial_lattice_vector_3: Optional[List[float]] = field(
-        default_factory=lambda: None
-    )
-    """The third lattice vector on the initial step, in Å."""
-    initial_lattice_vector_1_length: Optional[float] = None
-    """The length of the first lattice vector on the initial step, in Å."""
-    initial_lattice_vector_2_length: Optional[float] = None
-    """The length of the second lattice vector on the initial step, in Å."""
-    initial_lattice_vector_3_length: Optional[float] = None
-    """The length of the third lattice vector on the initial step, in Å."""
-    initial_angle_alpha: Optional[float] = None
-    """The angle between the second and third lattice vectors on the initial step, in degrees."""
-    initial_angle_beta: Optional[float] = None
-    """The angle between the first and third lattice vectors on the initial step, in degrees."""
-    initial_angle_gamma: Optional[float] = None
-    """The angle between the first and second lattice vectors on the initial step, in degrees."""
+    lattice_vector_1: Optional[List[float]] = field(default_factory=lambda: None)
+    """The first lattice vector, in Å."""
+    lattice_vector_2: Optional[List[float]] = field(default_factory=lambda: None)
+    """The second lattice vector, in Å."""
+    lattice_vector_3: Optional[List[float]] = field(default_factory=lambda: None)
+    """The third lattice vector, in Å."""
+    lattice_vector_1_length: Optional[float] = None
+    """The length of the first lattice vector, in Å."""
+    lattice_vector_2_length: Optional[float] = None
+    """The length of the second lattice vector, in Å."""
+    lattice_vector_3_length: Optional[float] = None
+    """The length of the third lattice vector, in Å."""
+    angle_alpha: Optional[float] = None
+    """The angle between the second and third lattice vectors, in degrees."""
+    angle_beta: Optional[float] = None
+    """The angle between the first and third lattice vectors, in degrees."""
+    angle_gamma: Optional[float] = None
+    """The angle between the first and second lattice vectors, in degrees."""
+
+    # stoichiometry data folded into the structure model
+    ion_types: Optional[List[str]] = field(default_factory=lambda: None)
+    """The distinct types of ions in the system."""
+    num_ion_types: Optional[List[int]] = field(default_factory=lambda: None)
+    """The number of ions of each type in the system."""
+    num_ion_types_primitive: Optional[List[int]] = field(default_factory=lambda: None)
+    """The number of ions of each type in the primitive cell."""
+    formula: Optional[str] = None
+    """The chemical formula of the system, in the format {element}{count if count > 1 else ''}, e.g. A3B2CD4."""
+    compound: Optional[str] = None
+    """The name of the compound, in the format {element1}-{element2}-..., e.g. A-B-C."""
 
 
 @dataclass
diff --git a/src/py4vasp/_raw/write.py b/src/py4vasp/_raw/write.py
index 1ce3d436..a3b89a71 100644
--- a/src/py4vasp/_raw/write.py
+++ b/src/py4vasp/_raw/write.py
@@ -25,6 +25,11 @@ def write(h5f, raw_data, *, selection=None):
 def _write_dataset(h5f, target, data, valid_indices=None):
     if check.is_none(data):
         return
+    if check.is_none(target):
+        # The chosen source defines no path for this field (e.g. the "ion"
+        # dielectric-function source has no q_point), so there is nothing to write.
+        # Skip it instead of passing the unset VaspData(None) target to h5py as a key.
+        return
     if isinstance(target, Link):
         write(h5f, data, selection=target.source)
     elif isinstance(target, Length) or target in h5f:
diff --git a/src/py4vasp/_util/database.py b/src/py4vasp/_util/database.py
index 0064c083..1699a34b 100644
--- a/src/py4vasp/_util/database.py
+++ b/src/py4vasp/_util/database.py
@@ -233,9 +233,24 @@ def get_formula_and_compound(
     tuple[str, str, list[str], list[int], list[int]]
         The chemical formula and compound name, followed by unique ion types, their
         counts in the conventional cell, and their counts in the primitive cell.
+
+    When the element names are unavailable but the counts are known, the counts (and
+    their primitive reduction) are still returned; the name-derived formula, compound,
+    and unique-type list are ``None`` because they cannot be built without names.
     """
-    if ion_types is None or number_ion_types is None:
+    if number_ion_types is None:
         return None, None, None, None, None
+    if ion_types is None:
+        # No element names: keep the counts in their raw order (they cannot be
+        # aggregated or sorted by element) and reduce them to the primitive cell.
+        simple_numbers = list(number_ion_types)
+        return (
+            None,
+            None,
+            None,
+            simple_numbers,
+            get_primitive_ion_numbers(simple_numbers),
+        )
 
     formula_dict = {}
     # first sort, then count up numbers in case any ion type is non-unique
diff --git a/tests/calculation/test_band.py b/tests/calculation/test_band.py
index 604d4ea6..c8a67f97 100644
--- a/tests/calculation/test_band.py
+++ b/tests/calculation/test_band.py
@@ -8,6 +8,7 @@
 import pytest
 
 from py4vasp import exception
+from py4vasp._calculation._dispersion import DispersionHandler
 from py4vasp._calculation.band import _OCCUPATION_CUTOFF, Band, BandHandler
 from py4vasp._calculation.kpoint import Kpoint
 from py4vasp._calculation.projector import Projector
@@ -696,6 +697,17 @@ def _check_to_database(_band):
         _band.ref, "fermi_energy_argument", _band.ref.fermi_energy
     )
 
+    # dispersion is folded into the band model
+    dispersion = DispersionHandler.from_data(
+        _band.ref.raw_data.dispersion
+    ).to_database()
+    assert database_data.eigenvalue_min == dispersion.eigenvalue_min
+    assert database_data.eigenvalue_max == dispersion.eigenvalue_max
+    assert database_data.eigenvalue_min_up == dispersion.eigenvalue_min_up
+    assert database_data.eigenvalue_max_up == dispersion.eigenvalue_max_up
+    assert database_data.eigenvalue_min_down == dispersion.eigenvalue_min_down
+    assert database_data.eigenvalue_max_down == dispersion.eigenvalue_max_down
+
     if getattr(_band.ref, "num_occupied_bands", None) is not None:
         assert database_data.num_occupied_bands == _band.ref.num_occupied_bands
     elif (
@@ -740,11 +752,12 @@ def test_to_database_noncollinear_projectors(noncollinear_projectors):
 
 
 def test_dispatcher_to_database_default(single_band):
-    """Dispatcher._to_database() must return {selection_name: handler_result}."""
+    """Dispatcher._to_database() returns {quantity: {selection: handler_result}}."""
     result = single_band._to_database()
     assert isinstance(result, dict)
     assert "band" in result
-    assert isinstance(result["band"], Band_DB)
+    assert isinstance(result["band"], dict)
+    assert isinstance(result["band"]["default"], Band_DB)
 
 
 def test_factory_methods(raw_data, check_factory_methods):
diff --git a/tests/calculation/test_bandgap.py b/tests/calculation/test_bandgap.py
index 1e2b5109..b0bff145 100644
--- a/tests/calculation/test_bandgap.py
+++ b/tests/calculation/test_bandgap.py
@@ -411,11 +411,12 @@ def test_to_database_spin_polarized(spin_polarized_handler, Assert):
 
 
 def test_dispatcher_to_database(bandgap):
-    """Dispatcher._to_database() must return {selection_name: handler_result}."""
+    """Dispatcher._to_database() returns {quantity: {selection: handler_result}}."""
     result = bandgap._to_database()
     assert isinstance(result, dict)
     assert "bandgap" in result
-    assert isinstance(result["bandgap"], Bandgap_DB)
+    assert isinstance(result["bandgap"], dict)
+    assert isinstance(result["bandgap"]["default"], Bandgap_DB)
 
 
 def test_dispatcher_to_dict_matches_read(raw_data, Assert):
diff --git a/tests/calculation/test_current_density.py b/tests/calculation/test_current_density.py
index 90f15273..d8bf5469 100644
--- a/tests/calculation/test_current_density.py
+++ b/tests/calculation/test_current_density.py
@@ -7,8 +7,9 @@
 import pytest
 
 from py4vasp import exception
-from py4vasp._calculation.current_density import CurrentDensity
+from py4vasp._calculation.current_density import CurrentDensity, CurrentDensityHandler
 from py4vasp._calculation.structure import Structure
+from py4vasp._raw.data_db import CurrentDensity_DB
 
 
 @pytest.fixture(params=("all", "x", "y", "z"))
@@ -193,6 +194,26 @@ def test_incorrect_selection_raises_error(raw_data):
         current_density.to_quiver("foo", a=0)
 
 
+def test_to_database(raw_data, Assert):
+    raw_current = raw_data.current_density("all")
+    handler = CurrentDensityHandler.from_data(raw_current)
+    db_data = handler.to_database()
+    assert isinstance(db_data, CurrentDensity_DB)
+    # magnitude |j| across all perturbations and grid points
+    magnitudes = np.concatenate(
+        [
+            np.linalg.norm(np.array(array), axis=0).ravel()
+            for array in raw_current.current_density
+        ]
+    )
+    assert db_data.magnitude_min == pytest.approx(float(np.min(magnitudes)))
+    assert db_data.magnitude_max == pytest.approx(float(np.max(magnitudes)))
+    assert db_data.magnitude_mean == pytest.approx(float(np.mean(magnitudes)))
+    # grid shape reported as [nx, ny, nz] from an array of shape (3, nz, ny, nx)
+    nz, ny, nx = np.array(raw_current.current_density[0]).shape[1:]
+    assert db_data.grid_shape == [nx, ny, nz]
+
+
 def test_factory_methods(raw_data, check_factory_methods):
     data = raw_data.current_density("x")
     check_factory_methods(CurrentDensity, data)
diff --git a/tests/calculation/test_dispatch.py b/tests/calculation/test_dispatch.py
index 16d6056e..2234809b 100644
--- a/tests/calculation/test_dispatch.py
+++ b/tests/calculation/test_dispatch.py
@@ -415,18 +415,18 @@ def test_default_selection_uses_quantity_name_as_key(self):
             result = merge_to_database(
                 source, "energy", _FakeHandler.from_data, _FakeHandler.read
             )
-        assert result == {"energy": {"value": 42}}
+        assert result == {"energy": {"default": {"value": 42}}}
 
-    def test_named_selection_appends_selection_to_quantity_name(self):
+    def test_named_selection_becomes_inner_key(self):
         raw = {"value": 10}
         source = DictSource({("band", "kpoints_opt"): raw})
         with _patch_sources(["kpoints_opt"]):
             result = merge_to_database(
                 source, "band", _FakeHandler.from_data, _FakeHandler.read
             )
-        assert result == {"band_kpoints_opt": {"value": 10}}
+        assert result == {"band": {"kpoints_opt": {"value": 10}}}
 
-    def test_multiple_selections_produce_distinct_keys(self):
+    def test_multiple_selections_produce_distinct_inner_keys(self):
         raw_a = {"value": 1}
         raw_b = {"value": 2}
         source = DictSource({("dos", "a"): raw_a, ("dos", "b"): raw_b})
@@ -434,7 +434,7 @@ def test_multiple_selections_produce_distinct_keys(self):
             result = merge_to_database(
                 source, "dos", _FakeHandler.from_data, _FakeHandler.read
             )
-        assert result == {"dos_a": {"value": 1}, "dos_b": {"value": 2}}
+        assert result == {"dos": {"a": {"value": 1}, "b": {"value": 2}}}
 
     def test_leading_underscore_stripped_from_quantity_name(self):
         raw = {"value": 7}
@@ -453,28 +453,28 @@ def test_leading_underscore_stripped_with_named_selection(self):
             result = merge_to_database(
                 source, "_CONTCAR", _FakeHandler.from_data, _FakeHandler.read
             )
-        assert result == {"CONTCAR_sel": {"value": 3}}
+        assert result == {"CONTCAR": {"sel": {"value": 3}}}
 
-    def test_default_key_has_no_default_suffix(self):
+    def test_default_selection_is_inner_key(self):
         raw = {"value": 5}
         source = DataSource(raw)
         with _patch_sources(["default"]):
             result = merge_to_database(
                 source, "force", _FakeHandler.from_data, _FakeHandler.read
             )
-        assert "force_default" not in result
-        assert "force" in result
+        assert set(result) == {"force"}
+        assert set(result["force"]) == {"default"}
 
     def test_duplicate_source_results_collapse_to_default(self):
         # an in-memory DataSource yields the same data for every source, so the
-        # non-default selections must collapse into the single default key
+        # non-default selections must collapse into the single default inner key
         raw = {"value": 99}
         source = DataSource(raw)
         with _patch_sources(["default", "final", "exciton"]):
             result = merge_to_database(
                 source, "structure", _FakeHandler.from_data, _FakeHandler.read
             )
-        assert result == {"structure": {"value": 99}}
+        assert result == {"structure": {"default": {"value": 99}}}
 
     def test_forwards_extra_kwargs_to_handler(self):
         raw = {"value": 4}
@@ -487,7 +487,7 @@ def test_forwards_extra_kwargs_to_handler(self):
                 _FakeHandler.read_with_args,
                 scale=3,
             )
-        assert result == {"energy": {"value": 12}}
+        assert result == {"energy": {"default": {"value": 12}}}
 
 
 @dataclasses.dataclass
@@ -600,14 +600,14 @@ def to_database(self):
             source, "quantity", _PartialHandler.from_data, _PartialHandler.to_database
         )
         assert "quantity" in result
-        assert isinstance(result["quantity"], _PartialDB)
+        assert isinstance(result["quantity"]["default"], _PartialDB)
 
     def test_non_empty_dict_result_included(self):
         source = DataSource({"value": 1})
         result = merge_to_database(
             source, "quantity", _FakeHandler.from_data, _FakeHandler.read
         )
-        assert result == {"quantity": {"value": 1}}
+        assert result == {"quantity": {"default": {"value": 1}}}
 
 
 class TestSubstituteRemainingSelection:
diff --git a/tests/calculation/test_optics.py b/tests/calculation/test_optics.py
index c9731e3f..167e1c1c 100644
--- a/tests/calculation/test_optics.py
+++ b/tests/calculation/test_optics.py
@@ -395,7 +395,8 @@ def test_to_database_keyed_by_optics(visible):
     assert isinstance(result, dict)
     # the DataSource ignores the selection, so all sources collapse to a single "optics"
     assert "optics" in result
-    assert isinstance(result["optics"], Optics_DB)
+    assert isinstance(result["optics"], dict)
+    assert isinstance(result["optics"]["default"], Optics_DB)
     # keys are derived from "optics", never from the underlying "dielectric_function"
     assert not any(key.startswith("dielectric_function") for key in result)
 
diff --git a/tests/calculation/test_pair_correlation.py b/tests/calculation/test_pair_correlation.py
index f58a8649..99629261 100644
--- a/tests/calculation/test_pair_correlation.py
+++ b/tests/calculation/test_pair_correlation.py
@@ -2,9 +2,10 @@
 # Licensed under the Apache License 2.0 (http://www.apache.org/licenses/LICENSE-2.0)
 from unittest.mock import patch
 
+import numpy as np
 import pytest
 
-from py4vasp import exception
+from py4vasp import exception, raw
 from py4vasp._calculation.pair_correlation import (
     PairCorrelation,
     PairCorrelationHandler,
@@ -109,6 +110,39 @@ def test_to_database(pair_correlation, raw_data):
     assert db_data.distance_max == float(pair_correlation.ref.distances[-1])
 
 
+def _pair_correlation_from_total(distances, total):
+    """Build a raw pair correlation whose only curve is the given total g(r)."""
+    function = np.asarray(total)[np.newaxis, np.newaxis, :]  # (steps, labels, points)
+    return raw.PairCorrelation(
+        distances=np.asarray(distances, dtype=float),
+        function=function,
+        labels=("total",),
+    )
+
+
+def test_to_database_first_peak():
+    """The first peak is the first local maximum of the total g(r) above the
+    threshold; a small sub-threshold bump before it must be ignored."""
+    distances = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
+    # bump at index 1 (0.5) is below the threshold; the first real peak is at index 4
+    total = [0.0, 0.5, 0.3, 2.5, 4.0, 1.5, 1.2, 1.0]
+    raw_pcf = _pair_correlation_from_total(distances, total)
+    db_data = PairCorrelationHandler.from_data(raw_pcf).to_database()
+    assert isinstance(db_data, PairCorrelation_DB)
+    assert db_data.first_peak_position == 4.0
+    assert db_data.first_peak_height == 4.0
+
+
+def test_to_database_no_first_peak():
+    """A monotonic curve staying below the threshold has no detectable peak."""
+    distances = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
+    total = [0.0, 0.2, 0.4, 0.6, 0.8, 0.9]
+    raw_pcf = _pair_correlation_from_total(distances, total)
+    db_data = PairCorrelationHandler.from_data(raw_pcf).to_database()
+    assert db_data.first_peak_position is None
+    assert db_data.first_peak_height is None
+
+
 def test_factory_methods(raw_data, check_factory_methods):
     data = raw_data.pair_correlation("Sr2TiO4")
     check_factory_methods(PairCorrelation, data)
diff --git a/tests/calculation/test_phonon_band.py b/tests/calculation/test_phonon_band.py
index 2fa52bb6..5e40dc31 100644
--- a/tests/calculation/test_phonon_band.py
+++ b/tests/calculation/test_phonon_band.py
@@ -6,9 +6,11 @@
 import numpy as np
 import pytest
 
+from py4vasp._calculation._dispersion import DispersionHandler
 from py4vasp._calculation._stoichiometry import Stoichiometry
 from py4vasp._calculation.kpoint import Kpoint
 from py4vasp._calculation.phonon_band import PhononBand, PhononBandHandler
+from py4vasp._raw.data_db import PhononBand_DB
 from py4vasp._util import convert
 
 
@@ -134,8 +136,14 @@ def test_print(phonon_band, format_):
 
 def test_to_database(phonon_band):
     handler = PhononBandHandler.from_data(phonon_band.ref.raw_data)
-    db_dict = handler.to_database()
-    assert db_dict == {}
+    db_data = handler.to_database()
+    assert isinstance(db_data, PhononBand_DB)
+    # dispersion (phonon frequencies) is folded into the phonon band model
+    dispersion = DispersionHandler.from_data(
+        phonon_band.ref.raw_data.dispersion
+    ).to_database()
+    assert db_data.eigenvalue_min == dispersion.eigenvalue_min
+    assert db_data.eigenvalue_max == dispersion.eigenvalue_max
 
 
 def test_factory_methods(raw_data, check_factory_methods):
diff --git a/tests/calculation/test_run_info.py b/tests/calculation/test_run_info.py
index 647d75e3..2c057a9f 100644
--- a/tests/calculation/test_run_info.py
+++ b/tests/calculation/test_run_info.py
@@ -118,11 +118,12 @@ def test_to_database(run_info_handler):
 
 
 def test_dispatcher_to_database(run_info):
-    """Dispatcher._to_database() must return {selection_name: handler_result}."""
+    """Dispatcher._to_database() returns {quantity: {selection: handler_result}}."""
     result = run_info._to_database()
     assert isinstance(result, dict)
     assert "run_info" in result
-    assert isinstance(result["run_info"], RunInfo_DB)
+    assert isinstance(result["run_info"], dict)
+    assert isinstance(result["run_info"]["default"], RunInfo_DB)
 
 
 def test_factory_methods(raw_data, check_factory_methods):
diff --git a/tests/calculation/test_stoichiometry.py b/tests/calculation/test_stoichiometry.py
index ae25c47d..70f3e956 100644
--- a/tests/calculation/test_stoichiometry.py
+++ b/tests/calculation/test_stoichiometry.py
@@ -218,9 +218,11 @@ def _setup(self, request, raw_data):
                 "text/html": "<em>A</em><sub>2</sub><em>B</em><em>C</em><sub>4</sub>",
             }
         self.ref_ion_types = ["O", "Sr", "Ti"] if request.param == "Sr2TiO4" else None
-        self.ref_num_ion_types = [4, 2, 1] if request.param == "Sr2TiO4" else None
+        # without element names the counts survive in their raw order (no name-based
+        # aggregation is possible), so num_ion_types is [2, 1, 4] rather than None
+        self.ref_num_ion_types = [4, 2, 1] if request.param == "Sr2TiO4" else [2, 1, 4]
         self.ref_num_ion_types_primitive = (
-            None if not request.param == "Sr2TiO4" else [4, 2, 1]
+            [4, 2, 1] if request.param == "Sr2TiO4" else [2, 1, 4]
         )
         self.ref_formula = None if not request.param == "Sr2TiO4" else "O4Sr2Ti"
         self.ref_compound = None if not request.param == "Sr2TiO4" else "O-Sr-Ti"
@@ -261,17 +263,15 @@ def test_factory_methods(raw_data, check_factory_methods):
     check_factory_methods(Stoichiometry, data, skip_methods=["to_mdtraj"])
 
 
-def test_to_database_dispatch(raw_data):
-    """The dispatcher keys the handler result by quantity name for the database."""
+def test_stoichiometry_not_collected_standalone(raw_data):
+    """Stoichiometry is folded into structure, so the dispatcher exposes no
+    ``_to_database`` and is skipped by the calculation-level collection. The
+    handler ``to_database`` stays (it is used by the structure model)."""
     raw_stoichiometry = raw_data.stoichiometry("Sr2TiO4")
-    result = Stoichiometry.from_data(raw_stoichiometry)._to_database()
-    assert set(result) == {"stoichiometry"}
-    expected = StoichiometryHandler.from_data(raw_stoichiometry).to_database()
-    assert result["stoichiometry"] == expected
-
-
-def test_to_database_dispatch_skips_empty(raw_data):
-    """Stoichiometry without ion types carries no database data, so it is omitted."""
-    raw_stoichiometry = raw_data.stoichiometry("Sr2TiO4 without ion types")
-    result = Stoichiometry.from_data(raw_stoichiometry)._to_database()
-    assert result == {}
+    dispatcher = Stoichiometry.from_data(raw_stoichiometry)
+    assert not hasattr(dispatcher, "_to_database")
+    # the handler still produces the stoichiometry model for the parent to fold in
+    assert isinstance(
+        StoichiometryHandler.from_data(raw_stoichiometry).to_database(),
+        Stoichiometry_DB,
+    )
diff --git a/tests/calculation/test_structure.py b/tests/calculation/test_structure.py
index 6942fa7e..1ba7e9bd 100644
--- a/tests/calculation/test_structure.py
+++ b/tests/calculation/test_structure.py
@@ -7,7 +7,7 @@
 import pytest
 
 from py4vasp import exception, raw
-from py4vasp._calculation._stoichiometry import Stoichiometry
+from py4vasp._calculation._stoichiometry import Stoichiometry, StoichiometryHandler
 from py4vasp._calculation.structure import Structure, StructureHandler
 from py4vasp._raw.data_db import Structure_DB
 from py4vasp._util import check
@@ -546,90 +546,112 @@ def test_system_dimensionality(Graphite, Sr2TiO4, Fe3O4):
 
 def test_to_database(structures, Assert):
     handler = StructureHandler.from_data(structures.ref.raw_data)
-    db_data: Structure_DB = handler.to_database()
-    assert isinstance(db_data, Structure_DB)
     has_timesteps = structures.ref.positions.ndim == 3
-    final_positions = (
+    dimensionality = structures.ref.dimensionality
+    num_ions = len(
         structures.ref.positions[-1] if has_timesteps else structures.ref.positions
     )
-    final_volume = structures.ref.volume[-1] if has_timesteps else structures.ref.volume
-    initial_volume = (
-        structures.ref.volume[0] if has_timesteps else structures.ref.volume
-    )
-    final_area_2d = (
-        (structures.ref.area_2d[-1] if structures.ref.area_2d is not None else None)
-        if has_timesteps
-        else structures.ref.area_2d
-    )
-    initial_area_2d = (
-        (structures.ref.area_2d[0] if structures.ref.area_2d is not None else None)
-        if has_timesteps
-        else structures.ref.area_2d
-    )
-    final_lattice_vectors = (
-        structures.ref.lattice_vectors[-1]
-        if has_timesteps
-        else structures.ref.lattice_vectors
-    )
-    initial_lattice_vectors = (
-        structures.ref.lattice_vectors[0]
-        if has_timesteps
-        else structures.ref.lattice_vectors
-    )
-    final_dimensionality = structures.ref.dimensionality
-    assert db_data.num_ions == len(final_positions)
-    assert db_data.dimensionality == final_dimensionality
-
-    for lattice_vectors, area_2d, volume, prefix in [
-        (final_lattice_vectors, final_area_2d, final_volume, "final"),
-        (initial_lattice_vectors, initial_area_2d, initial_volume, "initial"),
-    ]:
-        assert getattr(db_data, f"{prefix}_cell_volume") == pytest.approx(volume)
-        if final_dimensionality == 2:
-            assert getattr(db_data, f"{prefix}_cell_area_2d") == pytest.approx(area_2d)
-            assert getattr(db_data, f"{prefix}_cell_area_2d_span") == "12"
-        else:
-            assert getattr(db_data, f"{prefix}_cell_area_2d") is None
-            assert getattr(db_data, f"{prefix}_cell_area_2d_span") is None
 
-        for idx in range(3):
-            Assert.allclose(
-                getattr(db_data, f"{prefix}_lattice_vector_{idx+1}"),
-                lattice_vectors[idx],
+    def reference_geometry(index):
+        lattice = (
+            structures.ref.lattice_vectors[index]
+            if has_timesteps
+            else structures.ref.lattice_vectors
+        )
+        volume = (
+            structures.ref.volume[index] if has_timesteps else structures.ref.volume
+        )
+        area_2d = (
+            (
+                structures.ref.area_2d[index]
+                if structures.ref.area_2d is not None
+                else None
             )
+            if has_timesteps
+            else structures.ref.area_2d
+        )
+        return lattice, volume, area_2d
+
+    def check_geometry(db_data, index):
+        # Each model now holds a single geometry with unprefixed fields.
+        assert isinstance(db_data, Structure_DB)
+        assert db_data.num_ions == num_ions
+        assert db_data.dimensionality == dimensionality
+        lattice, volume, area_2d = reference_geometry(index)
+        assert db_data.cell_volume == pytest.approx(volume)
+        if dimensionality == 2:
+            assert db_data.cell_area_2d == pytest.approx(area_2d)
+            assert db_data.cell_area_2d_span == "12"
+        else:
+            assert db_data.cell_area_2d is None
+            assert db_data.cell_area_2d_span is None
+        for idx in range(3):
+            Assert.allclose(getattr(db_data, f"lattice_vector_{idx + 1}"), lattice[idx])
             assert getattr(
-                db_data, f"{prefix}_lattice_vector_{idx+1}_length"
-            ) == pytest.approx(np.linalg.norm(lattice_vectors[idx]))
-    alpha, beta, gamma = (
-        np.degrees(
-            np.arccos(
-                np.dot(lattice_vectors[i], lattice_vectors[j])
-                / (
-                    np.linalg.norm(lattice_vectors[i])
-                    * np.linalg.norm(lattice_vectors[j])
+                db_data, f"lattice_vector_{idx + 1}_length"
+            ) == pytest.approx(np.linalg.norm(lattice[idx]))
+        alpha, beta, gamma = (
+            np.degrees(
+                np.arccos(
+                    np.dot(lattice[i], lattice[j])
+                    / (np.linalg.norm(lattice[i]) * np.linalg.norm(lattice[j]))
                 )
             )
+            for i, j in ((1, 2), (0, 2), (0, 1))
         )
-        for i, j in ((1, 2), (0, 2), (0, 1))
-    )
-    assert getattr(db_data, f"{prefix}_angle_alpha") == pytest.approx(alpha)
-    assert getattr(db_data, f"{prefix}_angle_beta") == pytest.approx(beta)
-    assert getattr(db_data, f"{prefix}_angle_gamma") == pytest.approx(gamma)
+        assert db_data.angle_alpha == pytest.approx(alpha)
+        assert db_data.angle_beta == pytest.approx(beta)
+        assert db_data.angle_gamma == pytest.approx(gamma)
+
+    # the default describes the final geometry; steps=0 describes the initial one
+    check_geometry(handler.to_database(), -1)
+    check_geometry(handler.to_database(steps=0), 0)
+
+    # stoichiometry is folded into the structure model
+    stoichiometry = StoichiometryHandler.from_data(
+        structures.ref.raw_data.stoichiometry
+    ).to_database()
+    db_data = handler.to_database()
+    assert db_data.formula == stoichiometry.formula
+    assert db_data.compound == stoichiometry.compound
+    assert db_data.ion_types == stoichiometry.ion_types
+    assert db_data.num_ion_types == stoichiometry.num_ion_types
+    assert db_data.num_ion_types_primitive == stoichiometry.num_ion_types_primitive
+
+
+def test_to_database_split_distinct_initial_final(raw_data):
+    """A trajectory with distinct first/last steps yields both a final and an
+    initial structure model, each describing a single geometry."""
+    raw_structure = raw_data.structure("Fe3O4")  # trajectory with shifting positions
+    result = Structure.from_data(raw_structure)._to_database()
+    assert set(result) == {"structure"}
+    assert set(result["structure"]) == {"final", "initial"}
+    expected_final = StructureHandler.from_data(raw_structure).to_database(steps=-1)
+    expected_initial = StructureHandler.from_data(raw_structure).to_database(steps=0)
+    assert result["structure"]["final"] == expected_final
+    assert result["structure"]["initial"] == expected_initial
+    assert result["structure"]["final"] != result["structure"]["initial"]
+
+
+def test_to_database_identical_steps_only_final(raw_data):
+    """When the initial geometry matches the final one, only the final model is
+    stored (initial is dropped as a duplicate)."""
+    raw_structure = raw_data.structure("Sr2TiO4")  # trajectory of identical steps
+    result = Structure.from_data(raw_structure)._to_database()
+    assert set(result["structure"]) == {"final"}
 
 
-def test_to_database_dispatch(raw_data):
-    """The dispatcher keys the handler result by quantity name for the database."""
-    raw_structure = raw_data.structure("Sr2TiO4")
+def test_to_database_single_geometry_only_final(raw_data):
+    """A single-geometry structure is stored under the final key only."""
+    raw_structure = raw_data.structure("ZnS")  # not a trajectory
     result = Structure.from_data(raw_structure)._to_database()
-    assert set(result) == {"structure"}
-    expected = StructureHandler.from_data(raw_structure).to_database()
-    assert result["structure"] == expected
+    assert set(result["structure"]) == {"final"}
 
 
 def test_to_database_collects_available_sources(tmp_path):
     """Database collection enumerates every structure source. For a CONTCAR-only
-    calculation only the ``final`` source carries data, so the result is keyed
-    ``structure_final`` (default/exciton/poscar have no data and are dropped)."""
+    calculation only the ``final`` source carries data, so the result is nested as
+    ``{"structure": {"final": ...}}`` (default/exciton/poscar have no data)."""
     import h5py
 
     import py4vasp
@@ -639,8 +661,9 @@ def test_to_database_collects_available_sources(tmp_path):
         py4vasp._raw.write.write(h5f, raw.Version(99, 99, 99))
         py4vasp._raw.write.write(h5f, contcar_demo.Sr2TiO4())
     result = py4vasp.Calculation.from_path(tmp_path).structure._to_database()
-    assert set(result) == {"structure_final"}
-    assert result["structure_final"].num_ions == 7
+    assert set(result) == {"structure"}
+    assert set(result["structure"]) == {"final"}
+    assert result["structure"]["final"].num_ions == 7
 
 
 def test_factory_methods(raw_data, check_factory_methods):
diff --git a/tests/raw/test_write.py b/tests/raw/test_write.py
index 860be6de..a0adaa5f 100644
--- a/tests/raw/test_write.py
+++ b/tests/raw/test_write.py
@@ -45,6 +45,23 @@ def test_write_encodes_unicode_strings(tmp_path):
         assert np.array_equal(actual.number_ion_types, [2, 1, 4])
 
 
+def test_write_skips_field_absent_from_source_schema(tmp_path, Assert):
+    # The ionic() demo sets a q_point, but the "ion" source schema defines no q_point
+    # path. The writer must skip that field (this source simply does not store it)
+    # instead of crashing when it tries to use the unset schema target as an HDF5 key.
+    from py4vasp._demo.dielectric_function import ionic
+
+    filename = tmp_path / DEFAULT_FILE
+    raw_dielectric = ionic()
+    assert not raw_dielectric.q_point.is_none()  # the demo does set q_point
+    with h5py.File(filename, "w") as h5f:
+        write(h5f, raw.Version(99, 99, 99))
+        write(h5f, raw_dielectric, selection="ion")  # must not raise
+    with raw.access("dielectric_function", path=tmp_path, selection="ion") as actual:
+        Assert.allclose(actual.energies, raw_dielectric.energies)
+        assert actual.q_point.is_none()  # not stored by the "ion" source
+
+
 def test_write_skips_absent_mapping_entries(tmp_path):
     # A Mapping field is a per-index list whose entries may be absent (VaspData(None))
     # for the chosen selection. electron_phonon self_energy with "CRTA" has nbands_sum
diff --git a/tests/util/test_database.py b/tests/util/test_database.py
index b2d30212..be6a3651 100644
--- a/tests/util/test_database.py
+++ b/tests/util/test_database.py
@@ -130,6 +130,12 @@ def test_get_primitive_ion_numbers(ion_numbers, expected):
             [2, 3, 1, 1, 2],
             ("AsBrCa", "As-Br-Ca", ["As", "Br", "Ca"], [3, 3, 3], [1, 1, 1]),
         ),
+        # names absent but counts present: keep the counts (and primitive counts),
+        # but formula/compound/types cannot be derived without element names
+        (None, [2, 1, 4], (None, None, None, [2, 1, 4], [2, 1, 4])),
+        (None, [4, 2, 6], (None, None, None, [4, 2, 6], [2, 1, 3])),
+        # neither names nor counts: everything is None
+        (None, None, (None, None, None, None, None)),
     ],
 )
 def test_get_formula_and_compound(ion_types, ion_numbers, expectations):
@@ -189,7 +195,7 @@ def test_get_all_possible_keys():
     assert output_type_dict["band"] == output_type_dict["band:kpoints_opt"]
     assert output_type_dict["band"] == output_type_dict["band:kpoints_wan"]
     assert "band:default" not in output_type_dict
-    assert output_type_dict["current_density"] is None
+    assert output_type_dict["current_density"] == "CurrentDensity_DB"
 
     assert (
         sum([1 for v in all_keys.values() if len(v) > 0 and isinstance(v[0], tuple)])
@@ -207,9 +213,14 @@ def basic_db_checks(demo_calc_db: _DatabaseData, minimum_counter=1):
     # Check metadata fields
     assert isinstance(demo_calc_db.metadata.path, Path)
 
-    # Check that properties has enough non-empty entries
+    # properties is a dict of dicts {quantity: {selection: model}}; count the
+    # non-empty inner models across all quantities and selections.
     non_empty_counter = sum(
-        1 for v in demo_calc_db.properties.values() if v not in (None, {}, [])
+        1
+        for selections in demo_calc_db.properties.values()
+        if isinstance(selections, dict)
+        for v in selections.values()
+        if v not in (None, {}, [])
     )
     assert non_empty_counter > minimum_counter
     assert "run_info" in demo_calc_db.properties
diff --git a/tests/util/test_to_database_new.py b/tests/util/test_to_database_new.py
index 4a28be3d..ace496a0 100644
--- a/tests/util/test_to_database_new.py
+++ b/tests/util/test_to_database_new.py
@@ -7,8 +7,10 @@
   - CalculationMetaData has a `path` (directory) and `schema_version`,
     but no `tags` and no `hdf5_original_path`.
   - _to_database() takes no arguments.
-  - properties keys follow the format <quantity> (default) or
-    <quantity>_<selection> (non-default), with no leading underscore.
+  - properties is a dict of dicts: {"<quantity>": {"<selection>": model}}.
+    Top-level keys are bare quantity names (no leading underscore, no selection
+    suffix); the inner dict is keyed by selection with the default source keyed
+    "default".
   - schema_version is stored only in metadata, not on individual _DB dataclasses.
 """
 
@@ -119,6 +121,18 @@ def test_properties_has_entries(demo_db):
     assert len(demo_db.properties) > 0
 
 
+def test_properties_values_are_dicts(demo_db):
+    """properties is a dict of dicts: each value maps selection -> model."""
+    for key, value in demo_db.properties.items():
+        assert isinstance(value, dict), f"properties[{key!r}] is not a dict"
+        assert len(value) > 0, f"properties[{key!r}] is empty"
+
+
+def test_default_selection_is_inner_key(demo_db):
+    """run_info is always available under its default selection key 'default'."""
+    assert "default" in demo_db.properties["run_info"]
+
+
 def test_no_leading_underscore_in_properties_keys(demo_db):
     """Private quantities like _CONTCAR must be stored under 'CONTCAR', not '_CONTCAR'."""
     for key in demo_db.properties:
@@ -130,6 +144,27 @@ def test_run_info_in_properties(demo_db):
     assert "run_info" in demo_db.properties
 
 
+def test_subcomponents_not_top_level(demo_db):
+    """stoichiometry and dispersion are folded into their parent quantities and
+    must not appear as top-level properties."""
+    assert "stoichiometry" not in demo_db.properties
+    assert "dispersion" not in demo_db.properties
+
+
+def test_stoichiometry_folded_into_structure(demo_db):
+    """Structure models carry the folded stoichiometry fields."""
+    structure_models = list(demo_db.properties["structure"].values())
+    assert any(model.formula is not None for model in structure_models)
+    assert any(model.ion_types is not None for model in structure_models)
+
+
+def test_dispersion_folded_into_band(demo_db):
+    """The band model carries the folded dispersion eigenvalue range."""
+    band = demo_db.properties["band"]["default"]
+    assert band.eigenvalue_min is not None
+    assert band.eigenvalue_max is not None
+
+
 def test_default_selection_key_has_no_suffix(demo_db):
     """Keys for the default selection must not have a '_default' suffix."""
     for key in demo_db.properties:
@@ -137,14 +172,15 @@ def test_default_selection_key_has_no_suffix(demo_db):
 
 
 def test_non_default_selection_key_format(tmp_path):
-    """Non-default selections are appended with an underscore: quantity_selection."""
+    """Selections are nested keys, not folded into the top-level quantity key."""
     actual_path = tmp_path / "demo_calc_band"
     calc = demo.calculation(actual_path)
     db = calc._to_database()
-    # band has kpoints_opt and kpoints_wan selections in the schema.
-    # If the demo data doesn't have them, they simply won't appear — that is fine.
-    # But if they DO appear, the key format must be 'band_kpoints_opt' not
-    # 'band:kpoints_opt' or 'band.kpoints_opt'.
-    for key in db.properties:
-        assert ":" not in key, f"Key {key!r} contains a colon"
-        assert "." not in key, f"Key {key!r} contains a dot"
+    # Top-level keys are bare quantity names; selections live in the inner dict.
+    for quantity, selections in db.properties.items():
+        assert ":" not in quantity, f"Key {quantity!r} contains a colon"
+        assert "." not in quantity, f"Key {quantity!r} contains a dot"
+        assert "_" != quantity[-1:], f"Key {quantity!r} ends with underscore"
+        assert isinstance(selections, dict)
+        for selection in selections:
+            assert ":" not in selection, f"Selection {selection!r} contains a colon"