diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index de585175385c46080183b649f8d07a63798abf10..98d45eee9b9c3d4754573034bc0304943d662e95 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -179,7 +179,9 @@ class ParameterFile:
         "GEARFeedback": {
             "supernovae_energy_erg": 1e51,
             "yields_table": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
-            "discrete_yields": 0
+            "discrete_yields": 0,
+            "elements": ["Fe", "Mg", "O", "S", "Zn", "Sr", "Y", "Ba", "Eu"],
+            "metallicity_max_first_stars": -1,
         },
     }
     """
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index dbe718c9f17cf5c7f5ef01f313b5d694bd60ca0c..3c91a9f2c4e1b391d876f5da1be612d17ab48025 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -154,20 +154,19 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
 
       hydro_set_physical_internal_energy_dt(&p, &cosmo, 0);
       for(int j = 0; j < GEAR_CHEMISTRY_ELEMENT_COUNT; j++) {
-        p.chemistry_data.smoothed_metal_mass_fraction[j] = p.chemistry_data.metal_mass_fraction[j];
+        p.chemistry_data.smoothed_metal_mass_fraction[j] = p.chemistry_data.metal_mass[j] / hydro_get_mass(&p);
       }
 
 
       /* compute cooling rate */
       cooling_cool_part(&pconst, &us, &cosmo, &hydro_props,
-			&floor_props, &cooling, &p, &xp, /* Time */0, dt, dt);
+                        &floor_props, &cooling, &p, &xp, dt, dt, /* Time */ 0);
       float *tmp = PyArray_GETPTR1(rate, i);
       *tmp = hydro_get_physical_internal_energy_dt(&p, &cosmo);
-
     }
 
-  /* Cleanup */
-  cosmology_clean(&cosmo);
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
   cooling_clean(&cooling);
 
   return (PyObject *) rate;
diff --git a/src/pyswiftsim_tools.h b/src/pyswiftsim_tools.h
index 6b75075c21bc10bd6399cce584e128c83adb54ee..ac074fe8f977d74d84e1d898911a2c2b0c09a5f9 100644
--- a/src/pyswiftsim_tools.h
+++ b/src/pyswiftsim_tools.h
@@ -44,32 +44,39 @@ enum error_code {
 
 int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name);
 
-#define PYSWIFTSIM_INIT_STRUCTS()				\
-  								\
-  /* Initialize parser */					\
-  struct swift_params params;					\
-  parser_read_file(filename, &params);				\
-								\
-  /* Init unit system */					\
-  struct unit_system us;					\
+#define PYSWIFTSIM_INIT_STRUCTS()                             \
+                                                              \
+  /* Initialize parser */                                     \
+  struct swift_params params;                                 \
+  parser_read_file(filename, &params);                        \
+                                                              \
+  /* Init unit system */                                      \
+  struct unit_system us;                                      \
   units_init_from_params(&us, &params, "InternalUnitSystem");	\
-								\
-  /* Init physical constant */					\
-  struct phys_const pconst;					\
-  phys_const_init(&us, &params, &pconst);			\
-								\
-  /* Init cosmo */						\
-  struct cosmology cosmo;					\
-								\
-  if (with_cosmo)						\
-    cosmology_init(&params, &us, &pconst, &cosmo);		\
-  else								\
-    cosmology_init_no_cosmo(&cosmo);				\
-								\
-  /* Init hydro prop */						\
-  struct hydro_props hydro_props;				\
+                                                              \
+  /* Init physical constant */                            \
+  struct phys_const pconst;                               \
+  phys_const_init(&us, &params, &pconst);                 \
+                                                          \
+  /* Init cosmo */                                        \
+  struct cosmology cosmo;                                 \
+                                                          \
+  if (with_cosmo)                                         \
+    cosmology_init(&params, &us, &pconst, &cosmo);        \
+  else                                                    \
+    cosmology_init_no_cosmo(&cosmo);                      \
+                                                          \
+  /* Init hydro prop */                                   \
+  struct hydro_props hydro_props;                         \
   hydro_props_init(&hydro_props, &pconst, &us, &params);	\
-  
+
+
+
+#define PYSWIFTSIM_CLEAN_STRUCTS()              \
+  /* No need to clean swift_params */           \
+  /* No need to clean unit_system */            \
+  /* No need to clean phys_const */             \
+  cosmology_clean(&cosmo);                      \
 
 
 #endif // __PYSWIFTSIM_TOOLS_H__
diff --git a/src/stellar_evolution_wrapper.c b/src/stellar_evolution_wrapper.c
index 3a23c8d0fd80f6fbff27cc02c2b06e8fc267cf28..528e4c051f4e841922f78e664cb364b9f13b56fb 100644
--- a/src/stellar_evolution_wrapper.c
+++ b/src/stellar_evolution_wrapper.c
@@ -77,6 +77,10 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
     *imf_i /= pconst.const_solar_mass;
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) imf;
   
 }
@@ -142,6 +146,10 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
     *tau = pow(10, *tau) * pconst.const_year * 1e6;
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) lifetime;
   
 }
@@ -207,6 +215,10 @@ PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
     *m = pow(10, *m) * pconst.const_solar_mass;
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) mass;
   
 }
@@ -289,6 +301,10 @@ PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
     *r /= pconst.const_solar_mass * pconst.const_year * 1e6;
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) rate;
   
 }
@@ -370,6 +386,10 @@ PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
 
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) rate;
   
 }
@@ -443,6 +463,10 @@ PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args) {
     }
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) yields_out;
   
 }
@@ -510,6 +534,10 @@ PyObject* pysupernovae_ii_mass_ejected(PyObject* self, PyObject* args) {
 
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) mass_ejected;
   
 }
@@ -577,6 +605,10 @@ PyObject* pysupernovae_ii_mass_ejected_non_processed(PyObject* self, PyObject* a
     
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) mass_ejected;
   
 }
@@ -622,6 +654,10 @@ PyObject* pysupernovae_ia_yields(PyObject* self, PyObject* args) {
     *y = yields[i] / fp.stellar_model.snia.mass_white_dwarf;
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) yields_out;
   
 }
@@ -663,6 +699,10 @@ PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args) {
 
   *y = supernovae_ia_get_ejected_mass_processed(&fp.stellar_model.snia) / fp.stellar_model.snia.mass_white_dwarf;
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return (PyObject *) mass_ejected;
   
 }
@@ -707,6 +747,10 @@ PyObject* pystellar_get_elements(PyObject* self, PyObject* args) {
     }
   }
 
+  /* Clean the memory */
+  PYSWIFTSIM_CLEAN_STRUCTS();
+  feedback_clean(&fp);
+
   return t;
   
 }