diff --git a/nomad.yaml b/nomad.yaml
index bc68652c..8be3d8d5 100644
--- a/nomad.yaml
+++ b/nomad.yaml
@@ -6,7 +6,8 @@ plugins:
   # We only include our schema here. Without the explicit include, all plugins will be
   # loaded. Many build in plugins require more dependencies. Install nomad-lab[parsing]
   # to make all default plugins work.
-  include: 'schemas/nomad_simulations'
+  include:
+    - 'schemas/nomad_simulations'
   options:
     schemas/nomad_simulations:
       python_package: nomad_simulations
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 4e356c0d..acffd8aa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -18,7 +18,7 @@ authors = [
 license = { text = "Apache-2.0" }
 requires-python = ">=3.9"
 dependencies = [
-    "nomad-lab>=1.2.0",
+    "nomad-lab>=1.2.2dev578",
     "matid>=2.0.0.dev2",
 ]
 
diff --git a/src/nomad_simulations/outputs.py b/src/nomad_simulations/outputs.py
index d0ea0868..ff154a66 100644
--- a/src/nomad_simulations/outputs.py
+++ b/src/nomad_simulations/outputs.py
@@ -26,7 +26,9 @@
 from .model_system import ModelSystem
 from .physical_property import PhysicalProperty
 from .numerical_settings import SelfConsistency
-from .properties import ElectronicBandGap
+from .properties import (
+    ElectronicBandGap,
+)
 
 
 class Outputs(ArchiveSection):
@@ -62,8 +64,12 @@ class Outputs(ArchiveSection):
     # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
     # List of properties
     # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+
     electronic_band_gap = SubSection(sub_section=ElectronicBandGap.m_def, repeats=True)
 
+    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+
     def extract_spin_polarized_property(self, property_name: str) -> list:
         """
         Extracts the spin-polarized properties if present from the property name and returns them as a list of two elements in
@@ -84,9 +90,27 @@ def extract_spin_polarized_property(self, property_name: str) -> list:
             spin_polarized_properties.append(prop)
         return spin_polarized_properties
 
+    def set_model_system_ref(self) -> Optional[ModelSystem]:
+        """
+        Set the reference to the last ModelSystem if this is not set in the output. This is only
+        valid if there is only one ModelSystem in the parent section.
+
+        Returns:
+            (Optional[ModelSystem]): The reference to the last ModelSystem.
+        """
+        if self.m_parent is not None:
+            model_systems = self.m_parent.model_system
+            if model_systems is not None and len(model_systems) == 1:
+                return model_systems[-1]
+        return None
+
     def normalize(self, archive, logger) -> None:
         super().normalize(archive, logger)
 
+        # Set ref to the last ModelSystem if this is not set in the output
+        if self.model_system_ref is None:
+            self.model_system_ref = self.set_model_system_ref()
+
 
 class SCFOutputs(Outputs):
     """
diff --git a/src/nomad_simulations/physical_property.py b/src/nomad_simulations/physical_property.py
index f61a5e05..b013fa5b 100644
--- a/src/nomad_simulations/physical_property.py
+++ b/src/nomad_simulations/physical_property.py
@@ -171,7 +171,7 @@ def variables_shape(self) -> Optional[list]:
             (list): The shape of the variables over which the physical property varies.
         """
         if self.variables is not None:
-            return [v.get_n_grid_points(v.grid_points, logger) for v in self.variables]
+            return [v.get_n_grid_points(logger) for v in self.variables]
         return []
 
     @property
diff --git a/src/nomad_simulations/properties/band_gap.py b/src/nomad_simulations/properties/band_gap.py
index 87a9dd02..e1a424b7 100644
--- a/src/nomad_simulations/properties/band_gap.py
+++ b/src/nomad_simulations/properties/band_gap.py
@@ -33,7 +33,6 @@ class ElectronicBandGap(PhysicalProperty):
 
     iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandGap'
 
-    # ? can `type` change character depending on the `variables`?
     type = Quantity(
         type=MEnum('direct', 'indirect'),
         description="""
@@ -82,7 +81,7 @@ def __init__(
         self.name = self.m_def.name
         self.rank = []
 
-    def check_negative_values(self, logger: BoundLogger) -> Optional[pint.Quantity]:
+    def _check_negative_values(self, logger: BoundLogger) -> Optional[pint.Quantity]:
         """
         Checks if the electronic band gaps is negative and sets them to None if they are.
 
@@ -142,7 +141,7 @@ def normalize(self, archive, logger) -> None:
         super().normalize(archive, logger)
 
         # Checks if the `value` is negative and sets it to None if it is.
-        self.value = self.check_negative_values(logger)
+        self.value = self._check_negative_values(logger)
         if self.value is None:
             # ? What about deleting the class if `value` is None?
             logger.error('The `value` of the electronic band gap is not stored.')
diff --git a/src/nomad_simulations/variables.py b/src/nomad_simulations/variables.py
index 28c40a70..d5202d4a 100644
--- a/src/nomad_simulations/variables.py
+++ b/src/nomad_simulations/variables.py
@@ -55,38 +55,35 @@ class Variables(ArchiveSection):
 
     # grid_points_error = Quantity()
 
-    def get_n_grid_points(
-        self, grid_points: Optional[list], logger: BoundLogger
-    ) -> Optional[int]:
+    def get_n_grid_points(self, logger: BoundLogger) -> Optional[int]:
         """
         Get the number of grid points from the `grid_points` list. If `n_grid_points` is previously defined
         and does not coincide with the length of `grid_points`, a warning is issued and this function re-assigns `n_grid_points`
         as the length of `grid_points`.
 
         Args:
-            grid_points (Optional[list]): The grid points in which the variable is discretized.
             logger (BoundLogger): The logger to log messages.
 
         Returns:
             (Optional[int]): The number of grid points.
         """
-        if grid_points is not None and len(grid_points) > 0:
+        if self.grid_points is not None and len(self.grid_points) > 0:
             if (
-                self.n_grid_points != len(grid_points)
+                self.n_grid_points != len(self.grid_points)
                 and self.n_grid_points is not None
             ):
                 logger.warning(
                     f'The stored `n_grid_points`, {self.n_grid_points}, does not coincide with the length of `grid_points`, '
-                    f'{len(grid_points)}. We will re-assign `n_grid_points` as the length of `grid_points`.'
+                    f'{len(self.grid_points)}. We will re-assign `n_grid_points` as the length of `grid_points`.'
                 )
-            return len(grid_points)
+            return len(self.grid_points)
         return self.n_grid_points
 
     def normalize(self, archive, logger) -> None:
         super().normalize(archive, logger)
 
         # Setting `n_grid_points` if these are not defined
-        self.n_grid_points = self.get_n_grid_points(self.grid_points, logger)
+        self.n_grid_points = self.get_n_grid_points(logger)
 
 
 class Temperature(Variables):
diff --git a/tests/test_band_gap.py b/tests/test_band_gap.py
index 88ea1d47..c31f024b 100644
--- a/tests/test_band_gap.py
+++ b/tests/test_band_gap.py
@@ -58,7 +58,7 @@ def test_check_negative_values(
         self, value: Union[List[float], float], result: float
     ):
         """
-        Test the `check_negative_values` method.
+        Test the `_check_negative_values` method.
         """
         if isinstance(value, list):
             electronic_band_gap = ElectronicBandGap(
@@ -67,7 +67,7 @@ def test_check_negative_values(
         else:
             electronic_band_gap = ElectronicBandGap()
         electronic_band_gap.value = value * ureg.joule
-        checked_value = electronic_band_gap.check_negative_values(logger)
+        checked_value = electronic_band_gap._check_negative_values(logger)
         if checked_value is not None:
             assert np.isclose(checked_value.magnitude, result)
         else:
diff --git a/tests/test_variables.py b/tests/test_variables.py
index df4fbaab..f3bc4540 100644
--- a/tests/test_variables.py
+++ b/tests/test_variables.py
@@ -47,6 +47,6 @@ def test_normalize(self, n_grid_points: int, grid_points: list, result: int):
             n_grid_points=n_grid_points,
             grid_points=grid_points,
         )
-        assert variable.get_n_grid_points(grid_points, logger) == result
+        assert variable.get_n_grid_points(logger) == result
         variable.normalize(None, logger)
         assert variable.n_grid_points == result