diff --git a/festim/concentration/mobile.py b/festim/concentration/mobile.py
index f619e0f94..0af4900ee 100644
--- a/festim/concentration/mobile.py
+++ b/festim/concentration/mobile.py
@@ -61,7 +61,7 @@ def create_diffusion_form(
             raise ValueError(msg)
 
         F = 0
-        for material in materials.materials:
+        for material in materials:
             D_0 = material.D_0
             E_D = material.E_D
             c_0, c_0_n = self.get_concentration_for_a_given_material(material, T)
diff --git a/festim/concentration/theta.py b/festim/concentration/theta.py
index 11505eb72..51fb81615 100644
--- a/festim/concentration/theta.py
+++ b/festim/concentration/theta.py
@@ -30,7 +30,7 @@ def initialise(self, V, value, label=None, time_step=None):
         v = f.TestFunction(V)
         dx = f.Measure("dx", subdomain_data=self.volume_markers)
         F = 0
-        for mat in self.materials.materials:
+        for mat in self.materials:
             S = mat.S_0 * f.exp(-mat.E_S / k_B / self.T.T)
             F += -prev_sol * v * dx(mat.id)
             if mat.solubility_law == "sievert":
@@ -110,7 +110,7 @@ def create_form_post_processing(self, V, materials, dx):
         c = f.TrialFunction(V)
 
         F += -c * v * dx
-        for mat in materials.materials:
+        for mat in materials:
             if mat.solubility_law == "sievert":
                 # for sievert materials c = theta * S
                 F += self.solution * self.S * v * dx(mat.id)
diff --git a/festim/concentration/traps/trap.py b/festim/concentration/traps/trap.py
index d9916cabc..a1ad81b9c 100644
--- a/festim/concentration/traps/trap.py
+++ b/festim/concentration/traps/trap.py
@@ -151,7 +151,7 @@ def create_trapping_form(self, mobile, materials, T, dx, dt=None):
             # if the sim is steady state and
             # if a trap is not defined in one subdomain
             # add c_t = 0 to the form in this subdomain
-            for mat in materials.materials:
+            for mat in materials:
                 if mat not in self.materials:
                     F_trapping += solution * test_function * dx(mat.id)
 
diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py
index 03c4fd76a..8ca7c82bf 100644
--- a/festim/generic_simulation.py
+++ b/festim/generic_simulation.py
@@ -132,12 +132,12 @@ def materials(self):
 
     @materials.setter
     def materials(self, value):
-        if isinstance(value, list):
+        if isinstance(value, festim.Materials):
+            self._materials = value
+        elif isinstance(value, list):
             self._materials = festim.Materials(value)
         elif isinstance(value, festim.Material):
             self._materials = festim.Materials([value])
-        elif isinstance(value, festim.Materials):
-            self._materials = value
         elif value is None:
             self._materials = value
         else:
diff --git a/festim/materials/materials.py b/festim/materials/materials.py
index ea5c24fef..9d900b5b7 100644
--- a/festim/materials/materials.py
+++ b/festim/materials/materials.py
@@ -7,15 +7,13 @@
 from typing import Union
 
 
-class Materials:
+class Materials(list):
     """
-    Args:
-        materials (list, optional): contains festim.Material objects.
-            Defaults to [].
+    A list of festim.Material objects
     """
 
-    def __init__(self, materials=[]):
-        self.materials = materials
+    def __init__(self, *args):
+        super().__init__(*args)
         self.D = None
         self.S = None
         self.thermal_cond = None
@@ -38,7 +36,7 @@ def check_borders(self, size):
             bool -- True if everything's alright
         """
         all_borders = []
-        for m in self.materials:
+        for m in self:
             if isinstance(m.borders[0], list):
                 for border in m.borders:
                     all_borders.append(border)
@@ -63,7 +61,7 @@ def check_materials(self, T: festim.Temperature, derived_quantities: list = []):
                 objects the derived quantities. Defaults to [].
         """
 
-        if len(self.materials) > 0:  # TODO: get rid of this...
+        if len(self) > 0:  # TODO: get rid of this...
             self.check_consistency()
 
             self.check_for_unused_properties(T, derived_quantities)
@@ -75,7 +73,7 @@ def check_materials(self, T: festim.Temperature, derived_quantities: list = []):
     def check_unique_ids(self):
         # check that ids are different
         mat_ids = []
-        for mat in self.materials:
+        for mat in self:
             if type(mat.id) is list:
                 mat_ids += mat.id
             else:
@@ -98,12 +96,12 @@ def check_for_unused_properties(
         # warn about unused keys
         transient_properties = ["rho", "heat_capacity"]
         if not isinstance(T, HeatTransferProblem):
-            for mat in self.materials:
+            for mat in self:
                 for key in transient_properties:
                     if getattr(mat, key) is not None:
                         warnings.warn(key + " key will be ignored", UserWarning)
 
-        for mat in self.materials:
+        for mat in self:
             if getattr(mat, "thermal_cond") is not None:
                 warn = True
                 if isinstance(T, HeatTransferProblem):
@@ -135,9 +133,9 @@ def check_consistency(self):
         }
 
         for attr, value in attributes.items():
-            for mat in self.materials:
+            for mat in self:
                 value.append(getattr(mat, attr))
-            if value.count(None) not in [0, len(self.materials)]:
+            if value.count(None) not in [0, len(self)]:
                 raise ValueError("{} is not defined for all materials".format(attr))
 
     def check_missing_properties(self, T: festim.Temperature, derived_quantities: list):
@@ -151,12 +149,12 @@ def check_missing_properties(self, T: festim.Temperature, derived_quantities: li
             ValueError: if thermal_cond, heat_capacity or rho is None when needed
         """
         if isinstance(T, HeatTransferProblem):
-            if self.materials[0].thermal_cond is None:
+            if self[0].thermal_cond is None:
                 raise ValueError("Missing thermal_cond in materials")
             if T.transient:
-                if self.materials[0].heat_capacity is None:
+                if self[0].heat_capacity is None:
                     raise ValueError("Missing heat_capacity in materials")
-                if self.materials[0].rho is None:
+                if self[0].rho is None:
                     raise ValueError("Missing rho in materials")
         # TODO: add check for thermal cond for thermal flux computation
 
@@ -172,7 +170,7 @@ def find_material_from_id(self, mat_id):
         Returns:
             festim.Material: the material that has the id mat_id
         """
-        for material in self.materials:
+        for material in self:
             mat_ids = material.id
             if type(mat_ids) is not list:
                 mat_ids = [mat_ids]
@@ -192,7 +190,7 @@ def find_material_from_name(self, name):
         Returns:
             festim.Material: the material object
         """
-        for material in self.materials:
+        for material in self:
             if material.name == name:
                 return material
 
@@ -225,7 +223,7 @@ def find_subdomain_from_x_coordinate(self, x):
         Returns:
             int: the corresponding subdomain id
         """
-        for material in self.materials:
+        for material in self:
             # if no borders are provided, assume only one subdomain
             if material.borders is None:
                 return material.id
@@ -255,13 +253,13 @@ def create_properties(self, vm, T):
         """
         self.D = ArheniusCoeff(self, vm, T, "D_0", "E_D", degree=2)
         # all materials have the same properties so only checking the first is enough
-        if self.materials[0].S_0 is not None:
+        if self[0].S_0 is not None:
             self.S = ArheniusCoeff(self, vm, T, "S_0", "E_S", degree=2)
-        if self.materials[0].thermal_cond is not None:
+        if self[0].thermal_cond is not None:
             self.thermal_cond = ThermalProp(self, vm, T, "thermal_cond", degree=2)
             self.heat_capacity = ThermalProp(self, vm, T, "heat_capacity", degree=2)
             self.density = ThermalProp(self, vm, T, "rho", degree=2)
-        if self.materials[0].Q is not None:
+        if self[0].Q is not None:
             self.Q = ThermalProp(self, vm, T, "Q", degree=2)
 
     def solubility_as_function(self, mesh, T):
@@ -273,7 +271,7 @@ def solubility_as_function(self, mesh, T):
         vS = f.TestFunction(V)
         dx = mesh.dx
         F = 0
-        for mat in self.materials:
+        for mat in self:
             F += -S * vS * dx(mat.id)
             F += mat.S_0 * f.exp(-mat.E_S / k_B / T) * vS * dx(mat.id)
         f.solve(F == 0, S, bcs=[])
@@ -300,7 +298,7 @@ def create_solubility_law_markers(self, mesh: festim.Mesh):
         F_sievert = -sievert * test_function_sievert * mesh.dx
 
         # build the formulation depending on the
-        for mat in self.materials:
+        for mat in self:
             # make sure mat_ids is a list
             mat_ids = mat.id
             if not isinstance(mat.id, list):
diff --git a/festim/meshing/mesh_1d.py b/festim/meshing/mesh_1d.py
index ae8fb2219..f845b8972 100644
--- a/festim/meshing/mesh_1d.py
+++ b/festim/meshing/mesh_1d.py
@@ -73,7 +73,7 @@ def define_volume_markers(self, materials):
 
     def define_measures(self, materials):
         """Creates the fenics.Measure objects for self.dx and self.ds"""
-        if materials.materials[0].borders is not None:
+        if materials[0].borders is not None:
             materials.check_borders(self.size)
         self.define_markers(materials)
         super().define_measures()
diff --git a/festim/temperature/temperature_solver.py b/festim/temperature/temperature_solver.py
index b3111b398..263ba4529 100644
--- a/festim/temperature/temperature_solver.py
+++ b/festim/temperature/temperature_solver.py
@@ -122,7 +122,7 @@ def define_variational_problem(self, materials, mesh, dt=None):
         v_T = self.v_T
 
         self.F = 0
-        for mat in materials.materials:
+        for mat in materials:
             thermal_cond = mat.thermal_cond
             if callable(thermal_cond):  # if thermal_cond is a function
                 thermal_cond = thermal_cond(T)
diff --git a/test/system/test_misc.py b/test/system/test_misc.py
index e563b37e1..1dbc3fdbb 100644
--- a/test/system/test_misc.py
+++ b/test/system/test_misc.py
@@ -301,3 +301,14 @@ def test_small_timesteps_final_time_bug():
     my_model.run()
 
     assert np.isclose(my_model.t, my_model.settings.final_time, atol=0.0)
+
+
+def test_materials_setter():
+    """
+    Checks that @materials.setter properly assigns F.Materials to F.Simulation.materials
+    see #694 for the details
+    """
+    my_model = F.Simulation()
+    test_materials = F.Materials([])
+    my_model.materials = test_materials
+    assert my_model.materials is test_materials
diff --git a/test/unit/test_materials.py b/test/unit/test_materials.py
index c18d78b6d..3df541511 100644
--- a/test/unit/test_materials.py
+++ b/test/unit/test_materials.py
@@ -297,3 +297,13 @@ def test_error_wrong_solubility_law_string():
         match="Acceptable values for solubility_law are 'henry' and 'sievert'",
     ):
         F.Material(1, 1, 1, solubility_law="foo")
+
+
+def test_equality_identity_two_empty_materials():
+    """
+    Tests equality and two of two empty F.Materials objects, i.e. checks
+    that these F.Materials are equal but refer to different objects
+    """
+    my_materials1 = F.Materials([])
+    my_materials2 = F.Materials([])
+    assert (my_materials1 == my_materials2) and (my_materials1 is not my_materials2)