diff --git a/configure.ac b/configure.ac
index 2fa77df1c302d197eb83dc6f5b93a32181fcdb3d..c4a01a359ab83dca63cf29671639ecc0656c42d6 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1278,7 +1278,7 @@ case "$with_subgrid" in
 	with_subgrid_tracers=none
 	with_subgrid_entropy_floor=none
 	with_subgrid_stars=GEAR
-	with_subgrid_star_formation=none
+	with_subgrid_star_formation=GEAR
 	with_subgrid_feedback=thermal
    ;;
    EAGLE)
@@ -1787,7 +1787,7 @@ esac
 #  Star formation
 AC_ARG_WITH([star-formation], 
     [AS_HELP_STRING([--with-star-formation=<sfm>],
-       [star formation @<:@none, EAGLE, default: none@:>@] 
+       [star formation @<:@none, EAGLE, GEAR, default: none@:>@] 
     )],
     [with_star_formation="$withval"],
     [with_star_formation="none"]
@@ -1807,6 +1807,9 @@ case "$with_star_formation" in
    EAGLE)
       AC_DEFINE([STAR_FORMATION_EAGLE], [1], [EAGLE star formation model (Schaye and Dalla Vecchia (2008))])
    ;;
+   GEAR)
+      AC_DEFINE([STAR_FORMATION_GEAR], [1], [GEAR star formation model (Revaz and Jablonka (2018))])
+   ;;
    *)
       AC_MSG_ERROR([Unknown star formation model])
    ;;
diff --git a/src/Makefile.am b/src/Makefile.am
index 5a356ba6707eee0664b9cb7817c8f2195d5512b2..7ac5954ebf85818090065f50d6ab46ccb0786d62 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -146,6 +146,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 		 star_formation/none/star_formation_io.h \
 		 star_formation/EAGLE/star_formation.h star_formation/EAGLE/star_formation_struct.h \
 		 star_formation/EAGLE/star_formation_io.h \
+		 star_formation/GEAR/star_formation.h star_formation/GEAR/star_formation_struct.h \
+		 star_formation/GEAR/star_formation_io.h \
 		 cooling/none/cooling.h cooling/none/cooling_struct.h \
                  cooling/none/cooling_io.h \
 		 cooling/Compton/cooling.h cooling/Compton/cooling_struct.h \
diff --git a/src/runner.c b/src/runner.c
index cfe3398684bef7b5a04cbfcefce1ebfe7729e24b..83e489caba0140d941826fef3e8fedb094f212ea 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -62,6 +62,7 @@
 #include "space.h"
 #include "space_getsid.h"
 #include "star_formation.h"
+#include "star_formation_iact.h"
 #include "stars.h"
 #include "task.h"
 #include "timers.h"
@@ -1329,6 +1330,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
   const struct hydro_space *hs = &s->hs;
   const struct cosmology *cosmo = e->cosmology;
   const struct chemistry_global_data *chemistry = e->chemistry;
+  const struct star_formation *star_formation = e->star_formation;
   const float hydro_h_max = e->hydro_properties->h_max;
   const float hydro_h_min = e->hydro_properties->h_min;
   const float eps = e->hydro_properties->h_tolerance;
@@ -1412,6 +1414,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
           /* Finish the density calculation */
           hydro_end_density(p, cosmo);
           chemistry_end_density(p, chemistry, cosmo);
+          star_formation_end_density(p, star_formation, cosmo);
 
           /* Compute one step of the Newton-Raphson scheme */
           const float n_sum = p->density.wcount * h_old_dim;
@@ -1556,6 +1559,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             /* Re-initialise everything */
             hydro_init_part(p, hs);
             chemistry_init_part(p, chemistry);
+            star_formation_init_part(p, star_formation);
 
             /* Off we go ! */
             continue;
@@ -1574,6 +1578,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             if (has_no_neighbours) {
               hydro_part_has_no_neighbours(p, xp, cosmo);
               chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
+              star_formation_part_has_no_neighbours(p, xp, star_formation,
+                                                    cosmo);
             }
 
           } else {
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 383fa64c56f6e6e7a134452f9eef8f715171be0c..854f8b898df93be93bb020a46606a415170fe980 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -213,6 +213,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
       }
       if (r2 < hjg2 && pj_active) {
@@ -224,6 +225,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
         IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
       }
 
@@ -323,12 +325,14 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
           IACT(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
         } else if (pi_active) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
         } else if (pj_active) {
 
@@ -339,6 +343,7 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
         }
       }
@@ -426,12 +431,14 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
         IACT(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
       } else if (doi) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
       } else if (doj) {
 
@@ -442,6 +449,7 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
         IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -528,12 +536,14 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
         IACT(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
       } else if (doi) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
       } else if (doj) {
 
@@ -544,6 +554,7 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
         IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -629,6 +640,7 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
         IACT_NONSYM(r2, dx, hi, pj->h, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hi, pj->h, pi, pj, a, H);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -721,6 +733,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -776,6 +789,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -918,6 +932,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
         runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -1081,6 +1096,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -1168,6 +1184,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
         }
       } /* loop over the parts in ci. */
@@ -1457,6 +1474,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
         }
       } /* loop over the active parts in cj. */
@@ -1526,11 +1544,13 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
             IACT(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
           } else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
           }
         }
@@ -1625,6 +1645,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
         }
       } /* loop over the active parts in ci. */
@@ -1697,11 +1718,13 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
             IACT(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+            runner_iact_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
           } else {
             IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+            runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
           }
         }
@@ -1890,6 +1913,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
         }
       } /* loop over all other particles. */
@@ -1941,12 +1965,14 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
             IACT(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
           } else if (doi) {
 
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
           } else if (doj) {
 
@@ -1956,6 +1982,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
             IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+            runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
           }
         }
@@ -2082,6 +2109,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
           runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
 #endif
         }
       } /* loop over all other particles. */
@@ -2128,11 +2156,13 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
             IACT(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
           } else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
             runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
 #endif
           }
         }
diff --git a/src/space.c b/src/space.c
index 0d7621d4809662ec513c2717722d13871540470e..33e88a0b679b6b1500e5ce5ac657348843f46964 100644
--- a/src/space.c
+++ b/src/space.c
@@ -56,6 +56,7 @@
 #include "multipole.h"
 #include "restart.h"
 #include "sort_part.h"
+#include "star_formation.h"
 #include "stars.h"
 #include "threadpool.h"
 #include "tools.h"
@@ -3372,6 +3373,7 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
   const int with_gravity = e->policy & engine_policy_self_gravity;
 
   const struct chemistry_global_data *chemistry = e->chemistry;
+  const struct star_formation *star_formation = e->star_formation;
   const struct cooling_function_data *cool_func = e->cooling_func;
 
   /* Check that the smoothing lengths are non-zero */
@@ -3418,6 +3420,10 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
     /* Also initialise the chemistry */
     chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[k], &xp[k]);
 
+    /* Also initialise the star formation */
+    star_formation_first_init_part(phys_const, us, cosmo, star_formation, &p[k],
+                                   &xp[k]);
+
     /* And the cooling */
     cooling_first_init_part(phys_const, us, cosmo, cool_func, &p[k], &xp[k]);
 
diff --git a/src/star_formation.h b/src/star_formation.h
index 78b9e018bf5f6cabc8a2b34dc5c4915f12806f68..5f873d2da142ec6fda90e986b523f60f7ef0d4ef 100644
--- a/src/star_formation.h
+++ b/src/star_formation.h
@@ -32,6 +32,8 @@
 #include "./star_formation/none/star_formation.h"
 #elif defined(STAR_FORMATION_EAGLE)
 #include "./star_formation/EAGLE/star_formation.h"
+#elif defined(STAR_FORMATION_GEAR)
+#include "./star_formation/GEAR/star_formation.h"
 #else
 #error "Invalid choice of star formation law"
 #endif
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index e4753731b31df6634da9d3112369506fe4857ff5..b72bb38babaca51b3875147d04c46f2de95de1a7 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -626,4 +626,69 @@ INLINE static void starformation_print_backend(
           starform->max_gas_density_HpCM3);
 }
 
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Nothing to do here. We do not need to compute any quantity in the hydro
+ * density loop for the EAGLE star formation model.
+ *
+ * @param p The particle to act upon
+ * @param cd The global star_formation information.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_end_density(
+    struct part* restrict p, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * Nothing to do here. We do not need to compute any quantity in the hydro
+ * density loop for the EAGLE star formation model.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cd #star_formation containing star_formation informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_part_has_no_neighbours(struct part* restrict p,
+                                      struct xpart* restrict xp,
+                                      const struct star_formation* cd,
+                                      const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param cosmo The current cosmological model.
+ * @param data The global star_formation information used for this run.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_first_init_part(const struct phys_const* restrict phys_const,
+                               const struct unit_system* restrict us,
+                               const struct cosmology* restrict cosmo,
+                               const struct star_formation* data,
+                               const struct part* restrict p,
+                               struct xpart* restrict xp) {}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
+ *
+ * Nothing to do here. We do not need to compute any quantity in the hydro
+ * density loop for the EAGLE star formation model.
+ *
+ * @param p Pointer to the particle data.
+ * @param data The global star_formation information.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_init_part(
+    struct part* restrict p, const struct star_formation* data) {}
+
 #endif /* SWIFT_EAGLE_STAR_FORMATION_H */
diff --git a/src/star_formation/EAGLE/star_formation_iact.h b/src/star_formation/EAGLE/star_formation_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..84fbad9f230284d5679cc7dd304fd0dd376796a7
--- /dev/null
+++ b/src/star_formation/EAGLE/star_formation_iact.h
@@ -0,0 +1,69 @@
+/*******************************************************************************
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EAGLE_STAR_FORMATION_IACT_H
+#define SWIFT_EAGLE_STAR_FORMATION_IACT_H
+
+/**
+ * @file EAGLE/star_formation_iact.h
+ * @brief Density computation
+ */
+
+/**
+ * @brief do star_formation computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_star_formation(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
+
+  /* Nothing to do here. We do not need to compute any quantity in the hydro
+     density loop for the EAGLE star formation model. */
+}
+
+/**
+ * @brief do star_formation computation after the runner_iact_density (non
+ * symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj,
+                                  struct part *restrict pi,
+                                  const struct part *restrict pj, float a,
+                                  float H) {
+
+  /* Nothing to do here. We do not need to compute any quantity in the hydro
+     density loop for the EAGLE star formation model. */
+}
+
+#endif /* SWIFT_EAGLE_STAR_FORMATION_IACT_H */
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
new file mode 100644
index 0000000000000000000000000000000000000000..420cecb6f6cdfa6be8d0803064122fa539131d0c
--- /dev/null
+++ b/src/star_formation/GEAR/star_formation.h
@@ -0,0 +1,215 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *******************************************************************************/
+#ifndef SWIFT_GEAR_STAR_FORMATION_H
+#define SWIFT_GEAR_STAR_FORMATION_H
+
+/* Local includes */
+#include "cosmology.h"
+#include "error.h"
+#include "hydro_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Calculate if the gas has the potential of becoming
+ * a star.
+ *
+ * No star formation should occur, so return 0.
+ *
+ * @param starform the star formation law properties to use.
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param phys_const the physical constants in internal units.
+ * @param cosmo the cosmological parameters and properties.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cooling The cooling data struct.
+ *
+ */
+INLINE static int star_formation_is_star_forming(
+    const struct part* restrict p, const struct xpart* restrict xp,
+    const struct star_formation* starform, const struct phys_const* phys_const,
+    const struct cosmology* cosmo,
+    const struct hydro_props* restrict hydro_props,
+    const struct unit_system* restrict us,
+    const struct cooling_function_data* restrict cooling) {
+
+  return 0;
+}
+
+/**
+ * @brief Compute the star-formation rate of a given particle and store
+ * it into the #xpart.
+ *
+ * Nothing to do here.
+ *
+ * @param p #part.
+ * @param xp the #xpart.
+ * @param starform the star formation law properties to use
+ * @param phys_const the physical constants in internal units.
+ * @param cosmo the cosmological parameters and properties.
+ * @param dt_star The time-step of this particle.
+ */
+INLINE static void star_formation_compute_SFR(
+    const struct part* restrict p, struct xpart* restrict xp,
+    const struct star_formation* starform, const struct phys_const* phys_const,
+    const struct cosmology* cosmo, const double dt_star) {}
+
+/**
+ * @brief Decides whether a particle should be converted into a
+ * star or not.
+ *
+ * No SF should occur, so return 0.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param starform The properties of the star formation model.
+ * @param e The #engine (for random numbers).
+ * @param dt_star The time-step of this particle
+ * @return 1 if a conversion should be done, 0 otherwise.
+ */
+INLINE static int star_formation_should_convert_to_star(
+    const struct part* p, const struct xpart* xp,
+    const struct star_formation* starform, const struct engine* e,
+    const double dt_star) {
+
+  return 0;
+}
+
+/**
+ * @brief Update the SF properties of a particle that is not star forming.
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param e The #engine.
+ * @param starform The properties of the star formation model.
+ * @param with_cosmology Are we running with cosmology switched on?
+ */
+INLINE static void star_formation_update_part_not_SFR(
+    struct part* p, struct xpart* xp, const struct engine* e,
+    const struct star_formation* starform, const int with_cosmology) {}
+
+/**
+ * @brief Copies the properties of the gas particle over to the
+ * star particle.
+ *
+ * Nothing to do here.
+ *
+ * @param e The #engine
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param sp the new created star particle with its properties.
+ * @param starform the star formation law properties to use.
+ * @param phys_const the physical constants in internal units.
+ * @param cosmo the cosmological parameters and properties.
+ * @param with_cosmology if we run with cosmology.
+ */
+INLINE static void star_formation_copy_properties(
+    const struct part* p, const struct xpart* xp, struct spart* sp,
+    const struct engine* e, const struct star_formation* starform,
+    const struct cosmology* cosmo, const int with_cosmology) {}
+
+/**
+ * @brief initialization of the star formation law
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param starform the star formation law properties to initialize
+ *
+ */
+INLINE static void starformation_init_backend(
+    struct swift_params* parameter_file, const struct phys_const* phys_const,
+    const struct unit_system* us, const struct hydro_props* hydro_props,
+    const struct star_formation* starform) {}
+
+/**
+ * @brief Prints the used parameters of the star formation law
+ *
+ * @param starform the star formation law properties.
+ */
+INLINE static void starformation_print_backend(
+    const struct star_formation* starform) {
+
+  message("Star formation law is 'GEAR'");
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * @param p The particle to act upon
+ * @param cd The global star_formation information.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_end_density(
+    struct part* restrict p, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cd #star_formation containing star_formation informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_part_has_no_neighbours(struct part* restrict p,
+                                      struct xpart* restrict xp,
+                                      const struct star_formation* cd,
+                                      const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param cosmo The current cosmological model.
+ * @param data The global star_formation information used for this run.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_first_init_part(const struct phys_const* restrict phys_const,
+                               const struct unit_system* restrict us,
+                               const struct cosmology* restrict cosmo,
+                               const struct star_formation* data,
+                               const struct part* restrict p,
+                               struct xpart* restrict xp) {}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param data The global star_formation information.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_init_part(
+    struct part* restrict p, const struct star_formation* data) {}
+
+#endif /* SWIFT_GEAR_STAR_FORMATION_H */
diff --git a/src/star_formation/GEAR/star_formation_iact.h b/src/star_formation/GEAR/star_formation_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..ef4747001545e459380bed572ffb89848b649e5e
--- /dev/null
+++ b/src/star_formation/GEAR/star_formation_iact.h
@@ -0,0 +1,61 @@
+/*******************************************************************************
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GEAR_STAR_FORMATION_IACT_H
+#define SWIFT_GEAR_STAR_FORMATION_IACT_H
+
+/**
+ * @file GEAR/star_formation_iact.h
+ * @brief Density computation
+ */
+
+/**
+ * @brief do star_formation computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_star_formation(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {}
+
+/**
+ * @brief do star_formation computation after the runner_iact_density (non
+ * symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj,
+                                  struct part *restrict pi,
+                                  const struct part *restrict pj, float a,
+                                  float H) {}
+
+#endif /* SWIFT_GEAR_STAR_FORMATION_IACT_H */
diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ef04c49c4abcd00175aaa164271628a9ff89360
--- /dev/null
+++ b/src/star_formation/GEAR/star_formation_io.h
@@ -0,0 +1,44 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_STAR_FORMATION_GEAR_IO_H
+#define SWIFT_STAR_FORMATION_GEAR_IO_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "io_properties.h"
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param xparts The extended data particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+__attribute__((always_inline)) INLINE static int star_formation_write_particles(
+    const struct part* parts, const struct xpart* xparts,
+    struct io_props* list) {
+
+  return 0;
+}
+
+#endif /* SWIFT_STAR_FORMATION_GEAR_IO_H */
diff --git a/src/star_formation/GEAR/star_formation_struct.h b/src/star_formation/GEAR/star_formation_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..9b4e216fd2955f29d89dade6ee46c0e1af715cdb
--- /dev/null
+++ b/src/star_formation/GEAR/star_formation_struct.h
@@ -0,0 +1,31 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GEAR_STAR_FORMATION_STRUCT_H
+#define SWIFT_GEAR_STAR_FORMATION_STRUCT_H
+
+/**
+ * @brief Star-formation-related properties stored in the extended particle
+ * data.
+ */
+struct star_formation_xpart_data {};
+
+/* Starformation struct */
+struct star_formation {};
+
+#endif /* SWIFT_GEAR_STAR_FORMATION_STRUCT_H */
diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h
index 093c4118601ec79c4f55646975f0505b8357afa6..7dbe5b20cc401c0284036f3c973c9b65fcec8d2e 100644
--- a/src/star_formation/none/star_formation.h
+++ b/src/star_formation/none/star_formation.h
@@ -157,4 +157,62 @@ INLINE static void starformation_print_backend(
   message("Star formation law is 'No Star Formation'");
 }
 
+/**
+ * @brief Finishes the density calculation.
+ *
+ * @param p The particle to act upon
+ * @param cd The global star_formation information.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_end_density(
+    struct part* restrict p, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cd #star_formation containing star_formation informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_part_has_no_neighbours(struct part* restrict p,
+                                      struct xpart* restrict xp,
+                                      const struct star_formation* cd,
+                                      const struct cosmology* cosmo) {}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param cosmo The current cosmological model.
+ * @param data The global star_formation information used for this run.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_first_init_part(const struct phys_const* restrict phys_const,
+                               const struct unit_system* restrict us,
+                               const struct cosmology* restrict cosmo,
+                               const struct star_formation* data,
+                               const struct part* restrict p,
+                               struct xpart* restrict xp) {}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param data The global star_formation information.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_init_part(
+    struct part* restrict p, const struct star_formation* data) {}
+
 #endif /* SWIFT_NONE_STAR_FORMATION_H */
diff --git a/src/star_formation/none/star_formation_iact.h b/src/star_formation/none/star_formation_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd74115bec699029748806b512c9d6bd7fb829fe
--- /dev/null
+++ b/src/star_formation/none/star_formation_iact.h
@@ -0,0 +1,61 @@
+/*******************************************************************************
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_NONE_STAR_FORMATION_IACT_H
+#define SWIFT_NONE_STAR_FORMATION_IACT_H
+
+/**
+ * @file none/star_formation_iact.h
+ * @brief Density computation
+ */
+
+/**
+ * @brief do star_formation computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_star_formation(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {}
+
+/**
+ * @brief do star_formation computation after the runner_iact_density (non
+ * symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj,
+                                  struct part *restrict pi,
+                                  const struct part *restrict pj, float a,
+                                  float H) {}
+
+#endif /* SWIFT_NONE_STAR_FORMATION_IACT_H */
diff --git a/src/star_formation_iact.h b/src/star_formation_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..ef457214a23102bc33385705db41c89dc29d8b8f
--- /dev/null
+++ b/src/star_formation_iact.h
@@ -0,0 +1,41 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_STAR_FORMATION_IACT_H
+#define SWIFT_STAR_FORMATION_IACT_H
+
+/**
+ * @file src/star_formation_iact.h
+ * @brief Branches between the different star formation iact.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Import the right star formation law definition */
+#if defined(STAR_FORMATION_NONE)
+#include "./star_formation/none/star_formation_iact.h"
+#elif defined(STAR_FORMATION_EAGLE)
+#include "./star_formation/EAGLE/star_formation_iact.h"
+#elif defined(STAR_FORMATION_GEAR)
+#include "./star_formation/GEAR/star_formation_iact.h"
+#else
+#error "Invalid choice of star formation law"
+#endif
+
+#endif /* SWIFT_STAR_FORMATION_IACT_H */
diff --git a/src/star_formation_io.h b/src/star_formation_io.h
index 248866c513e405b80f974d4932f02d3c707785ee..6c57bdc819fb8b583d217dbd60f402918d8f81ef 100644
--- a/src/star_formation_io.h
+++ b/src/star_formation_io.h
@@ -32,6 +32,8 @@
 #include "./star_formation/none/star_formation_io.h"
 #elif defined(STAR_FORMATION_EAGLE)
 #include "./star_formation/EAGLE/star_formation_io.h"
+#elif defined(STAR_FORMATION_GEAR)
+#include "./star_formation/GEAR/star_formation_io.h"
 #else
 #error "Invalid choice of star formation model."
 #endif
diff --git a/src/star_formation_struct.h b/src/star_formation_struct.h
index 3e9859cfb650ddd370c6acce282084b7eba6bf82..2a62d284b435c353525311979b343754856364e8 100644
--- a/src/star_formation_struct.h
+++ b/src/star_formation_struct.h
@@ -32,6 +32,8 @@
 #include "./star_formation/none/star_formation_struct.h"
 #elif defined(STAR_FORMATION_EAGLE)
 #include "./star_formation/EAGLE/star_formation_struct.h"
+#elif defined(STAR_FORMATION_GEAR)
+#include "./star_formation/GEAR/star_formation_struct.h"
 #else
 #error "Invalid choice of star formation structure."
 #endif
diff --git a/src/tools.c b/src/tools.c
index 45b04e35bb4713a5b73158839e642d0796aa94a7..43ac0177daef171850ea325f9fa23770fb82ae13 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -45,6 +45,7 @@
 #include "part.h"
 #include "periodic.h"
 #include "runner.h"
+#include "star_formation_iact.h"
 #include "stars.h"
 
 /**
@@ -222,6 +223,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
         /* Interact */
         runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj, a, H);
         runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hi, pj->h, pi, pj, a, H);
       }
     }
   }
@@ -254,6 +256,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
         /* Interact */
         runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi, a, H);
         runner_iact_nonsym_chemistry(r2, dx, hj, pi->h, pj, pi, a, H);
+        runner_iact_nonsym_star_formation(r2, dx, hj, pi->h, pj, pi, a, H);
       }
     }
   }
@@ -522,6 +525,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
         /* Interact */
         runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj, a, H);
         runner_iact_nonsym_chemistry(r2, dxi, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_star_formation(r2, dxi, hi, hj, pi, pj, a, H);
       }
 
       /* Hit or miss? */
@@ -534,6 +538,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
         /* Interact */
         runner_iact_nonsym_density(r2, dxi, hj, hi, pj, pi, a, H);
         runner_iact_nonsym_chemistry(r2, dxi, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_star_formation(r2, dxi, hj, hi, pj, pi, a, H);
       }
     }
   }