diff --git a/ThermoScreening/thermo/thermo.py b/ThermoScreening/thermo/thermo.py
index 135d751..a6b13b2 100644
--- a/ThermoScreening/thermo/thermo.py
+++ b/ThermoScreening/thermo/thermo.py
@@ -132,6 +132,9 @@ def run(self):
         -------
         None
         """
+        # Disabled: _transform_units is the only method that mutates the shared
+        # System in place (frequencies/energy). temperature_scan relies on run()
+        # leaving the System untouched, so keep this off (or make it copy first).
         # self._transform_units()
         self._rotational_contribution()
         self._vibrational_contribution()
@@ -142,6 +145,41 @@ def run(self):
         self._summary()
 
 
+    def temperature_scan(self, temperatures):
+        """
+        Recompute the thermochemistry at each of ``temperatures`` (in Kelvin).
+
+        The electronic energy, geometry, and vibrational frequencies are
+        temperature-independent, so this reuses the same :class:`System` (and
+        this object's pressure, engine, and quasi-RRHO setting) and only the
+        thermal terms are recomputed -- i.e. a temperature scan from a single
+        Hessian. This object is left unchanged.
+
+        Parameters
+        ----------
+        temperatures : iterable of float
+            Temperatures in Kelvin.
+
+        Returns
+        -------
+        list of Thermo
+            A ``Thermo`` computed (``run()`` already called) for each
+            temperature, in the given order.
+        """
+        scan = []
+        for temperature in temperatures:
+            thermo = Thermo(
+                temperature=temperature,
+                pressure=self._pressure,
+                system=self._system,
+                engine=self._engine,
+                quasi_rrho=self._quasi_rrho,
+            )
+            thermo.run()
+            scan.append(thermo)
+        return scan
+
+
     def _transform_units(self):
         """
         Transforms the units of the system to the correct units.
diff --git a/docs/api.rst b/docs/api.rst
index 5259b64..0e187f4 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -43,7 +43,8 @@ Thermochemistry core
 
 .. autoclass:: ThermoScreening.thermo.thermo.Thermo
    :members: total_energy, total_enthalpy, total_gibbs_free_energy,
-             total_entropy, total_heat_capacity, total_EeGtot, electronic_energy
+             total_entropy, total_heat_capacity, total_EeGtot, electronic_energy,
+             temperature_scan
 
 Coordinate and frequency readers
 --------------------------------
diff --git a/docs/usage.rst b/docs/usage.rst
index 28a4def..da0b037 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -98,6 +98,22 @@ The ensemble free energy lies at or below the lowest conformer's by the mixing
 (conformational) entropy; use ``lowest_gibbs`` when you instead want the single
 dominant structure to carry forward (e.g. into :func:`reaction_free_energy`).
 
+Temperature scans
+-----------------
+
+The electronic energy, geometry and frequencies do not depend on temperature,
+so a single (expensive) DFTB+/Hessian run can be reused across many
+temperatures:
+
+.. code-block:: python
+
+   thermo = dftbplus_thermo(molecule("H2O"), **dftb_3ob_parameters)
+   scan = thermo.temperature_scan([250, 273.15, 298.15, 350, 400])
+   gibbs = [(t._temperature, t.total_gibbs_free_energy("H")) for t in scan]
+
+Each entry is a fully-evaluated ``Thermo`` at that temperature; the original
+object is left unchanged.
+
 Reactions and redox
 -------------------
 
diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py
index 19879c7..6954304 100644
--- a/tests/thermo/test_thermo.py
+++ b/tests/thermo/test_thermo.py
@@ -123,6 +123,30 @@ def test_thermo_rejects_negative_pressure():
         Thermo(temperature=298.15, pressure=-1.0, system=_valid_system(), engine="dftb+")
 
 
+def test_temperature_scan_reuses_system_and_varies_with_temperature():
+    system = _valid_system()
+    thermo = Thermo(temperature=298.15, pressure=101325, system=system, engine="dftb+")
+
+    temperatures = [250.0, 298.15, 350.0]
+    scan = thermo.temperature_scan(temperatures)
+
+    assert [t._temperature for t in scan] == temperatures
+    # the expensive System (geometry/frequencies) is reused, not rebuilt
+    assert all(t._system is system for t in scan)
+    # temperature actually changes the thermochemistry
+    gibbs = [t.total_gibbs_free_energy("H") for t in scan]
+    assert len(set(gibbs)) == 3
+    # the scan leaves this object untouched
+    assert thermo._temperature == 298.15
+
+
+def test_temperature_scan_empty_and_invalid():
+    thermo = Thermo(temperature=298.15, pressure=101325, system=_valid_system(), engine="dftb+")
+    assert thermo.temperature_scan([]) == []
+    with pytest.raises(TSValueError, match="temperature is negative"):
+        thermo.temperature_scan([300.0, -5.0])
+
+
 # --- geometry-dependent thermochemistry, validated against ASE IdealGasThermo --- #
 
 from ase import Atoms  # noqa: E402