diff --git a/src/drift.h b/src/drift.h
new file mode 100644
index 0000000000000000000000000000000000000000..73e909e31093f2461b78e7d378ba79aa0bc7e255
--- /dev/null
+++ b/src/drift.h
@@ -0,0 +1,106 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_DRIFT_H
+#define SWIFT_DRIFT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "debug.h"
+
+/**
+ * @brief Perform the 'drift' operation on a #gpart
+ *
+ * @param gp The #gpart to drift.
+ * @param dt The drift time-step
+ * @param timeBase The minimal allowed time-step size.
+ * @param ti_old Integer start of time-step
+ * @param ti_current Integer end of time-step
+ */
+__attribute__((always_inline)) INLINE static void drift_gpart(struct gpart* gp,
+ float dt,
+ double timeBase,
+ int ti_old,
+ int ti_current) {
+ /* Drift... */
+ gp->x[0] += gp->v_full[0] * dt;
+ gp->x[1] += gp->v_full[1] * dt;
+ gp->x[2] += gp->v_full[2] * dt;
+
+ /* Compute offset since last cell construction */
+ gp->x_diff[0] -= gp->v_full[0] * dt;
+ gp->x_diff[1] -= gp->v_full[1] * dt;
+ gp->x_diff[2] -= gp->v_full[2] * dt;
+}
+
+/**
+ * @brief Perform the 'drift' operation on a #part
+ *
+ * @param p The #part to drift.
+ * @param xp The #xpart of the particle.
+ * @param dt The drift time-step
+ * @param timeBase The minimal allowed time-step size.
+ * @param ti_old Integer start of time-step
+ * @param ti_current Integer end of time-step
+ */
+__attribute__((always_inline)) INLINE static void drift_part(struct part* p,
+ struct xpart* xp,
+ float dt,
+ double timeBase,
+ int ti_old,
+ int ti_current) {
+ /* Useful quantity */
+ const float h_inv = 1.0f / p->h;
+
+ /* Drift... */
+ p->x[0] += xp->v_full[0] * dt;
+ p->x[1] += xp->v_full[1] * dt;
+ p->x[2] += xp->v_full[2] * dt;
+
+ /* Predict velocities (for hydro terms) */
+ p->v[0] += p->a_hydro[0] * dt;
+ p->v[1] += p->a_hydro[1] * dt;
+ p->v[2] += p->a_hydro[2] * dt;
+
+ /* Predict smoothing length */
+ const float w1 = p->h_dt * h_inv * dt;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ /* Predict density */
+ const float w2 = -3.0f * p->h_dt * h_inv * dt;
+ if (fabsf(w2) < 0.2f)
+ p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+ else
+ p->rho *= expf(w2);
+
+ /* Predict the values of the extra fields */
+ hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
+
+ /* Compute offset since last cell construction */
+ xp->x_diff[0] -= xp->v_full[0] * dt;
+ xp->x_diff[1] -= xp->v_full[1] * dt;
+ xp->x_diff[2] -= xp->v_full[2] * dt;
+}
+
+#endif /* SWIFT_DRIFT_H */
diff --git a/src/runner.c b/src/runner.c
index e0a52517d22fc991e3891a7868bc76c1f43e431f..c382da9b8eaed017ed84636f17d04c09f7cf3c5a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -47,6 +47,7 @@
#include "gravity.h"
#include "hydro_properties.h"
#include "hydro.h"
+#include "drift.h"
#include "kick.h"
#include "minmax.h"
#include "scheduler.h"
@@ -757,14 +758,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* Drift... */
- gp->x[0] += gp->v_full[0] * dt;
- gp->x[1] += gp->v_full[1] * dt;
- gp->x[2] += gp->v_full[2] * dt;
-
- /* Compute offset since last cell construction */
- gp->x_diff[0] -= gp->v_full[0] * dt;
- gp->x_diff[1] -= gp->v_full[1] * dt;
- gp->x_diff[2] -= gp->v_full[2] * dt;
+ drift_gpart(gp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
@@ -781,40 +775,8 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
- /* Useful quantity */
- const float h_inv = 1.0f / p->h;
-
/* Drift... */
- p->x[0] += xp->v_full[0] * dt;
- p->x[1] += xp->v_full[1] * dt;
- p->x[2] += xp->v_full[2] * dt;
-
- /* Predict velocities (for hydro terms) */
- p->v[0] += p->a_hydro[0] * dt;
- p->v[1] += p->a_hydro[1] * dt;
- p->v[2] += p->a_hydro[2] * dt;
-
- /* Predict smoothing length */
- const float w1 = p->h_dt * h_inv * dt;
- if (fabsf(w1) < 0.2f)
- p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
- else
- p->h *= expf(w1);
-
- /* Predict density */
- const float w2 = -3.0f * p->h_dt * h_inv * dt;
- if (fabsf(w2) < 0.2f)
- p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
- else
- p->rho *= expf(w2);
-
- /* Predict the values of the extra fields */
- hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
-
- /* Compute offset since last cell construction */
- xp->x_diff[0] -= xp->v_full[0] * dt;
- xp->x_diff[1] -= xp->v_full[1] * dt;
- xp->x_diff[2] -= xp->v_full[2] * dt;
+ drift_part(p, xp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +