diff --git a/atomistics/calculators/ase.py b/atomistics/calculators/ase.py
index 4f58eb81..b446fc19 100644
--- a/atomistics/calculators/ase.py
+++ b/atomistics/calculators/ase.py
@@ -396,7 +396,7 @@ def calc_molecular_dynamics_thermal_expansion_with_ase(
     temperature_lst = np.arange(
         temperature_start, temperature_stop + temperature_step, temperature_step
     ).tolist()
-    volume_md_lst, temperature_md_lst = [], []
+    volume_md_lst, temperature_md_lst, cell_md_lst = [], [], []
     for temperature in get_tqdm_iterator(temperature_lst):
         result_dict = calc_molecular_dynamics_npt_with_ase(
             structure=structure_current.copy(),
@@ -410,11 +410,13 @@ def calc_molecular_dynamics_thermal_expansion_with_ase(
             externalstress=externalstress,
         )
         structure_current.set_cell(cell=result_dict["cell"][-1], scale_atoms=True)
+        cell_md_lst.append(result_dict["cell"][-1])
         temperature_md_lst.append(result_dict["temperature"][-1])
         volume_md_lst.append(result_dict["volume"][-1])
     return get_thermal_expansion_output(
-        temperatures_lst=temperature_md_lst,
-        volumes_lst=volume_md_lst,
+        cell_lst=np.array(cell_md_lst),
+        temperatures_lst=np.array(temperature_md_lst),
+        volumes_lst=np.array(volume_md_lst),
         output_keys=output_keys,
     )
 
diff --git a/atomistics/calculators/lammps/helpers.py b/atomistics/calculators/lammps/helpers.py
index 29eb4e22..f41b9831 100644
--- a/atomistics/calculators/lammps/helpers.py
+++ b/atomistics/calculators/lammps/helpers.py
@@ -120,7 +120,7 @@ def lammps_thermal_expansion_loop(
         **kwargs,
     )
 
-    volume_md_lst, temperature_md_lst = [], []
+    volume_md_lst, cell_md_lst, temperature_md_lst = [], [], []
     for temp in get_tqdm_iterator(temperature_lst):
         run_str_rendered = Template(run_str).render(
             run=run,
@@ -133,12 +133,14 @@ def lammps_thermal_expansion_loop(
         )
         for line in run_str_rendered.split("\n"):
             lmp_instance.interactive_lib_command(line)
+        cell_md_lst.append(lmp_instance.interactive_cells_getter())
         volume_md_lst.append(lmp_instance.interactive_volume_getter())
         temperature_md_lst.append(lmp_instance.interactive_temperatures_getter())
     lammps_shutdown(lmp_instance=lmp_instance, close_instance=lmp is None)
     return get_thermal_expansion_output(
-        temperatures_lst=temperature_md_lst,
-        volumes_lst=volume_md_lst,
+        cell_lst=np.array(cell_md_lst),
+        temperatures_lst=np.array(temperature_md_lst),
+        volumes_lst=np.array(volume_md_lst),
         output_keys=output_keys,
     )
 
diff --git a/atomistics/shared/output.py b/atomistics/shared/output.py
index 620e4de5..19bba0b0 100644
--- a/atomistics/shared/output.py
+++ b/atomistics/shared/output.py
@@ -55,6 +55,7 @@ class OutputMolecularDynamics(Output):
 class OutputThermalExpansion(Output):
     temperatures: callable
     volumes: callable
+    cells: callable
 
 
 @dataclasses.dataclass
diff --git a/atomistics/shared/thermal_expansion.py b/atomistics/shared/thermal_expansion.py
index b714c14a..23b0cbe9 100644
--- a/atomistics/shared/thermal_expansion.py
+++ b/atomistics/shared/thermal_expansion.py
@@ -4,7 +4,12 @@
 
 
 class ThermalExpansionProperties:
-    def __init__(self, temperatures_lst: np.ndarray, volumes_lst: np.ndarray):
+    def __init__(
+        self,
+        cell_lst: np.ndarray,
+        temperatures_lst: np.ndarray,
+        volumes_lst: np.ndarray,
+    ):
         """
         Initialize the ThermalExpansionProperties class.
 
@@ -12,9 +17,19 @@ def __init__(self, temperatures_lst: np.ndarray, volumes_lst: np.ndarray):
         temperatures_lst (np.ndarray): Array of temperatures.
         volumes_lst (np.ndarray): Array of volumes.
         """
+        self._cell_lst = cell_lst
         self._temperatures_lst = temperatures_lst
         self._volumes_lst = volumes_lst
 
+    def cells(self) -> np.ndarray:
+        """
+        Get the array of cells.
+
+        Returns:
+        np.ndarray: Array of cells.
+        """
+        return self._cell_lst
+
     def volumes(self) -> np.ndarray:
         """
         Get the array of volumes.
@@ -35,13 +50,17 @@ def temperatures(self) -> np.ndarray:
 
 
 def get_thermal_expansion_output(
-    temperatures_lst: np.ndarray, volumes_lst: np.ndarray, output_keys: tuple[str]
+    temperatures_lst: np.ndarray,
+    volumes_lst: np.ndarray,
+    cell_lst: np.ndarray,
+    output_keys: tuple[str],
 ) -> dict:
     """
     Get the thermal expansion output.
 
     Parameters:
     temperatures_lst (np.ndarray): Array of temperatures.
+    cells_lst (np.ndarray): Array of cells.
     volumes_lst (np.ndarray): Array of volumes.
     output_keys (tuple[str]): Tuple of output keys.
 
@@ -49,7 +68,7 @@ def get_thermal_expansion_output(
     dict: Dictionary containing the thermal expansion output.
     """
     thermal = ThermalExpansionProperties(
-        temperatures_lst=temperatures_lst, volumes_lst=volumes_lst
+        cell_lst=cell_lst, temperatures_lst=temperatures_lst, volumes_lst=volumes_lst
     )
     return OutputThermalExpansion(
         **{k: getattr(thermal, k) for k in OutputThermalExpansion.keys()}
diff --git a/atomistics/workflows/phonons/helper.py b/atomistics/workflows/phonons/helper.py
index 33c66864..20fc99e0 100644
--- a/atomistics/workflows/phonons/helper.py
+++ b/atomistics/workflows/phonons/helper.py
@@ -224,6 +224,18 @@ def volumes(self) -> np.ndarray:
             * len(self._thermal_properties["temperatures"])
         )
 
+    def cells(self) -> np.ndarray:
+        """
+        Get the cells.
+
+        Returns:
+            np.ndarray: The cells.
+        """
+        return np.array(
+            [self._phonopy.unitcell.cell]
+            * len(self._thermal_properties["temperatures"])
+        )
+
 
 def restore_magmoms(
     structure_with_magmoms: Atoms,