diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 1c96323e755a9b7c19292ce66c24e4a5943b626a..4458e4012c5ea343295a97903214f732f5b8a5a7 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -45,34 +45,18 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(struct gpar
 
 
 
-__attribute__((always_inline)) INLINE static float external_gravity_compute_timestep(
-    struct gpart* g) {
-
-  /* Currently no limit is imposed */
-  const float dx   = g->x[0]-External_Potential_X;
-  const float dy   = g->x[1]-External_Potential_Y;
-  const float dz   = g->x[2]-External_Potential_Z;
-  const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
-  const float drdv = (g->x[0]-External_Potential_X) * (g->v_full[0]) + (g->x[1]-External_Potential_Y) * (g->v_full[1]) + (g->x[2]-External_Potential_Z) * (g->v_full[2]); 
-  const float dota_x = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
-  const float dota_y = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
-  const float dota_z = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
-  const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
-  const float a_2    = g->a_grav_external[0] * g->a_grav_external[0] + g->a_grav_external[1] * g->a_grav_external[1] + g->a_grav_external[2] * g->a_grav_external[2];
-  
-  return 0.03f * sqrtf(a_2/dota_2);
+__attribute__((always_inline)) INLINE static float external_gravity_compute_timestep(struct gpart* g) {
+  float dtmin = FLT_MAX;
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+  dtmin = fminf(dtmin, external_gravity_pointmass_timestep(g));
+#endif
+  return dtmin;
 }
 
-__attribute__((always_inline)) INLINE static void external_gravity_pointmass(struct gpart *g, float* agrav_x, float* agrav_y, float* agrav_z)
+__attribute__((always_inline)) INLINE static void external_gravity(struct gpart *g)
 {
-  const float dx   = g->x[0]-External_Potential_X;
-  const float dy   = g->x[1]-External_Potential_Y;
-  const float dz   = g->x[2]-External_Potential_Z;
-  const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
-
-
-  *agrav_x = - const_G *  External_Potential_Mass * dx * rinv * rinv * rinv;
-  *agrav_y = - const_G *  External_Potential_Mass * dy * rinv * rinv * rinv;
-  *agrav_z = - const_G *  External_Potential_Mass * dz * rinv * rinv * rinv;
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+  external_gravity_pointmass(g);
+#endif
 
 }
diff --git a/src/gravity/Default/potentials.h b/src/gravity/Default/potentials.h
index 3c71dec4c583dd06e55a7fa01339669c93d47c16..39cde2d1dd50631615ab5474b15d85b4302b0746 100644
--- a/src/gravity/Default/potentials.h
+++ b/src/gravity/Default/potentials.h
@@ -11,4 +11,39 @@
 #define External_Potential_Mass (1e10 * SOLAR_MASS_IN_CGS / const_unit_mass_in_cgs)
 #endif
 
+
+__attribute__((always_inline)) INLINE static float external_gravity_pointmass_timestep(
+    struct gpart* g) {
+
+  /* Currently no limit is imposed */
+  const float dx   = g->x[0]-External_Potential_X;
+  const float dy   = g->x[1]-External_Potential_Y;
+  const float dz   = g->x[2]-External_Potential_Z;
+  const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
+  const float drdv = (g->x[0]-External_Potential_X) * (g->v_full[0]) + (g->x[1]-External_Potential_Y) * (g->v_full[1]) + (g->x[2]-External_Potential_Z) * (g->v_full[2]); 
+  const float dota_x = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
+  const float dota_y = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
+  const float dota_z = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
+  const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+  const float a_2    = g->a_grav_external[0] * g->a_grav_external[0] + g->a_grav_external[1] * g->a_grav_external[1] + g->a_grav_external[2] * g->a_grav_external[2];
+  
+  return 0.03f * sqrtf(a_2/dota_2);
+}
+
+__attribute__((always_inline)) INLINE static void external_gravity_pointmass(struct gpart *g)
+{
+  const float dx   = g->x[0]-External_Potential_X;
+  const float dy   = g->x[1]-External_Potential_Y;
+  const float dz   = g->x[2]-External_Potential_Z;
+  const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
+  
+
+  g->a_grav_external[0] += - const_G *  External_Potential_Mass * dx * rinv * rinv * rinv;
+  g->a_grav_external[1] += - const_G *  External_Potential_Mass * dy * rinv * rinv * rinv;
+  g->a_grav_external[2] += - const_G *  External_Potential_Mass * dz * rinv * rinv * rinv;
+}
+
+
+
 #endif /* SWIFT_GRAVITY_CONST_H */
+
diff --git a/src/runner.c b/src/runner.c
index 6f9307810ec55fbc615d885e9ed6229e73cf3b15..cd4ed426f3318bc5aa981d9fa07ca4d17bf1ebfc 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -129,13 +129,8 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
 
     /* Is this part within the time step? */
     if (g->ti_end <= ti_current) {
-		float agrav_x, agrav_y, agrav_z;
-		external_gravity_pointmass(g, &agrav_x, &agrav_y, &agrav_z);
-		g->a_grav_external[0] += agrav_x;
-		g->a_grav_external[1] += agrav_y;
-		g->a_grav_external[2] += agrav_z;
+		external_gravity_pointmass(g);
 		  
-
       /* check for energy and angular momentum conservation - begin by synchronizing velocity*/
       const float dx   = g->x[0]-External_Potential_X;
       const float dy   = g->x[1]-External_Potential_Y;