diff --git a/configure.ac b/configure.ac
index 2597d61d01b9df1fa22c1759960179c25b42e226..32ef711f2e306cf16edd6ac5f58a08bd2c172dc1 100644
--- a/configure.ac
+++ b/configure.ac
@@ -917,7 +917,26 @@ case "$with_subgrid" in
    ;;
 esac
 
+# Gravity scheme.
+AC_ARG_WITH([gravity],
+   [AS_HELP_STRING([--with-gravity=<scheme>],
+      [Gravity scheme to use @<:@default, with-potential, default: default@:>@]
+   )],
+   [with_gravity="$withval"],
+   [with_gravity="default"]
+)
 
+case "$with_gravity" in
+   with-potential)
+      AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
+   ;;
+   default)
+      AC_DEFINE([DEFAULT_GRAVITY], [1], [Default gravity scheme])
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown gravity scheme: $with_gravity])
+   ;;
+esac
 
 # Hydro scheme.
 AC_ARG_WITH([hydro],
@@ -1351,12 +1370,14 @@ AC_MSG_RESULT([
    Equation of state  : $with_eos
    Adiabatic index    : $with_gamma
    Riemann solver     : $with_riemann
-   Cooling function   : $with_cooling
-   Chemistry          : $with_chemistry
 
-   External potential  : $with_potential
+   Gravity scheme      : $with_gravity
    Multipole order     : $with_multipole_order
    No gravity below ID : $no_gravity_below_id
+   External potential  : $with_potential
+
+   Cooling function   : $with_cooling
+   Chemistry          : $with_chemistry
 
    Individual timers     : $enable_timers
    Task debugging        : $enable_task_debugging
diff --git a/examples/AgoraDisk/agora_disk.yml b/examples/AgoraDisk/agora_disk.yml
index 598f1db2858c3e420c55ae31c9ca10d37c2f48b7..7368700d8a2a5ca8de7d677e1da78be51d669835 100644
--- a/examples/AgoraDisk/agora_disk.yml
+++ b/examples/AgoraDisk/agora_disk.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085e21   # kpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
diff --git a/examples/AgoraDisk/plotSolution.py b/examples/AgoraDisk/plotSolution.py
index 2c3de1728eb6f22bee72361db1e7f0c77613eb25..fd5ed8a5e9fb1a1475202633595699d682edf281 100644
--- a/examples/AgoraDisk/plotSolution.py
+++ b/examples/AgoraDisk/plotSolution.py
@@ -82,8 +82,8 @@ mean_dispersion_time_range       = [50, 500]     # in Myr      = Used if add_mea
 mean_dispersion_density_range    = [1e-25, 1e-22]# in g/cm3    = Used if add_mean_fractional_dispersion = 1, for draw_density_DF
 disk_normal_vector               = [0., 0., 1.]
 
-gadget_default_unit_base = {'UnitLength_in_cm'         : 3.08568e+21,
-                            'UnitMass_in_g'            :   1.989e+43,
+gadget_default_unit_base = {'UnitLength_in_cm'         : 3.08567758e+21,
+                            'UnitMass_in_g'            :    1.98848e+43,
                             'UnitVelocity_in_cm_per_s' :      100000}
 color_names              = ['red', 'magenta', 'orange', 'gold', 'green', 'cyan', 'blue', 'blueviolet', 'black']
 linestyle_names          = ['-']
diff --git a/examples/DiscPatch/GravityOnly/disc-patch.yml b/examples/DiscPatch/GravityOnly/disc-patch.yml
index c76e4f612250d180f2ba2fccd0c6209878173433..4ec061add978bec82c267660cc343cf0bfa8f4c6 100644
--- a/examples/DiscPatch/GravityOnly/disc-patch.yml
+++ b/examples/DiscPatch/GravityOnly/disc-patch.yml
@@ -1,8 +1,8 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e33     # Grams
-  UnitLength_in_cgs:   3.0856776e18  # Centimeters
-  UnitVelocity_in_cgs: 1e5           # Centimeters per second
+  UnitMass_in_cgs:     1.98848e33    # M_sun in grams
+  UnitLength_in_cgs:   3.08567758e18 # parsec in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in cm/s
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
diff --git a/examples/DiscPatch/GravityOnly/makeIC.py b/examples/DiscPatch/GravityOnly/makeIC.py
index 6e4e148392eb7ca2fbf8c29c3f737d029916c59b..5f9650f44277cf858021c9b628d68134c47a19b7 100644
--- a/examples/DiscPatch/GravityOnly/makeIC.py
+++ b/examples/DiscPatch/GravityOnly/makeIC.py
@@ -38,11 +38,10 @@ import random
 # usage: python makeIC.py 1000 
 
 # physical constants in cgs
-NEWTON_GRAVITY_CGS  = 6.672e-8
-SOLAR_MASS_IN_CGS   = 1.9885e33
-PARSEC_IN_CGS       = 3.0856776e18
-PROTON_MASS_IN_CGS  = 1.6726231e24
-YEAR_IN_CGS         = 3.154e+7
+NEWTON_GRAVITY_CGS  = 6.67408e-8
+SOLAR_MASS_IN_CGS   = 1.98848e33
+PARSEC_IN_CGS       = 3.08567758e18
+YEAR_IN_CGS         = 3.15569252e7
 
 # choice of units
 const_unit_length_in_cgs   =   (PARSEC_IN_CGS)
diff --git a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
index 2859599c03c44ff6a15522164a8e2a407bdde41e..983a7dcc103135ab4db61d6ea77701532226c101 100644
--- a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
+++ b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
@@ -1,8 +1,8 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e33         # Grams
-  UnitLength_in_cgs:   3.08567758149e18  # Centimeters
-  UnitVelocity_in_cgs: 1e5               # Centimeters per second
+  UnitMass_in_cgs:     1.98848e33     # Grams
+  UnitLength_in_cgs:   3.08567758e18  # Centimeters
+  UnitVelocity_in_cgs: 1e5            # Centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
diff --git a/examples/DiscPatch/HydroStatic/disc-patch.yml b/examples/DiscPatch/HydroStatic/disc-patch.yml
index 8816bc17ca526d01b7abcf55bb43287bbb36224a..422e1cf910202e8f6dc0a9395fc7e36ce80443ed 100644
--- a/examples/DiscPatch/HydroStatic/disc-patch.yml
+++ b/examples/DiscPatch/HydroStatic/disc-patch.yml
@@ -1,8 +1,8 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e33         # Grams
-  UnitLength_in_cgs:   3.08567758149e18  # Centimeters
-  UnitVelocity_in_cgs: 1e5               # Centimeters per second
+  UnitMass_in_cgs:     1.98848e33     # Grams
+  UnitLength_in_cgs:   3.08567758e18  # Centimeters
+  UnitVelocity_in_cgs: 1e5            # Centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
index 83b07e50d59869c51fd3854a7f15cfd7d476d04d..8b4c55560c34e7bdb538f2b4732369216f91a087 100644
--- a/examples/DiscPatch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -54,8 +54,8 @@ fileName = "Disc-Patch.hdf5"
 
 # physical constants in cgs
 NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.9885e33
-PARSEC_IN_CGS       = 3.08567758149e18
+SOLAR_MASS_IN_CGS   = 1.98848e33
+PARSEC_IN_CGS       = 3.08567758e18
 PROTON_MASS_IN_CGS  = 1.672621898e-24
 BOLTZMANN_IN_CGS    = 1.38064852e-16
 YEAR_IN_CGS         = 3.15569252e7
diff --git a/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml b/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml
index 6f17cfbb1e0125faf8e47fe4e9e55bfdf4df7b71..450689034f4ae782cc74bf01dac93e723e5d2ce2 100644
--- a/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml
+++ b/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e33         # Grams
-  UnitLength_in_cgs:   3.08567758149e18  # Centimeters
+  UnitMass_in_cgs:     1.98848e33        # Grams
+  UnitLength_in_cgs:   3.08567758e18     # Centimeters
   UnitVelocity_in_cgs: 1e5               # Centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
diff --git a/examples/DiscPatch/HydroStatic_1D/makeIC.py b/examples/DiscPatch/HydroStatic_1D/makeIC.py
index 1589dfc8c73e5b9bf3c2cad4bcf3029654d9e67e..983a550a3442c6470611792081a5884d38023a6a 100644
--- a/examples/DiscPatch/HydroStatic_1D/makeIC.py
+++ b/examples/DiscPatch/HydroStatic_1D/makeIC.py
@@ -55,8 +55,8 @@ fileName = "Disc-Patch.hdf5"
 
 # physical constants in cgs
 NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.9885e33
-PARSEC_IN_CGS       = 3.08567758149e18
+SOLAR_MASS_IN_CGS   = 1.98848e33
+PARSEC_IN_CGS       = 3.08567758e18
 PROTON_MASS_IN_CGS  = 1.672621898e-24
 BOLTZMANN_IN_CGS    = 1.38064852e-16
 YEAR_IN_CGS         = 3.15569252e7
diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml
index 3183f40fb1a727167dad32b2c01f995fea0ae61b..439bb7eb6dc5d460752771addc83c89e27f69b7f 100644
--- a/examples/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_100/eagle_100.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
@@ -28,10 +28,9 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  scale_factor_first:  0.92  # Scale-factor of the first snaphot (cosmological run)
+  scale_factor_first:  0.91  # Scale-factor of the first snaphot (cosmological run)
   time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.10  # Time difference between consecutive outputs (in internal units)
-  compression:         1
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
diff --git a/examples/EAGLE_100/run.sh b/examples/EAGLE_100/run.sh
index 9310fc370acb51256b46e6ba0fb739fb2728caf9..9c990a902a6350eff348aad40c482723d1ba954c 100755
--- a/examples/EAGLE_100/run.sh
+++ b/examples/EAGLE_100/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -t 16 eagle_100.yml 2>&1 | tee output.log
+../swift -c -s -G -S -t 16 eagle_100.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index 321c62cdf7dca35381541b2c6230260da4b0dbe2..8ebe29fb0216e16aeaafcdc086085d8c9879fc5f 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
@@ -28,9 +28,9 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  scale_factor_first:  0.92  # Scale-factor of the first snaphot (cosmological run)
+  scale_factor_first:  0.91  # Scale-factor of the first snaphot (cosmological run)
   time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.10  # Time difference between consecutive outputs (in internal units)
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
   compression:         1
 
 # Parameters governing the conserved quantities statistics
diff --git a/examples/EAGLE_12/run.sh b/examples/EAGLE_12/run.sh
index 82d0e3e4d4418f0d00851f1279587a1ed45144df..67f1c24a1ead927823b9240cdeb718b35580d573 100755
--- a/examples/EAGLE_12/run.sh
+++ b/examples/EAGLE_12/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -t 16 eagle_12.yml 2>&1 | tee output.log
+../swift -c -s -G -S -t 16 eagle_12.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 682a5485e8319b9946aac5234075e1f7ba395a48..904f8c58168e0bae0ca36c5c1a71f1d006d4e3f2 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
@@ -28,9 +28,9 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  scale_factor_first:  0.92  # Scale-factor of the first snaphot (cosmological run)
+  scale_factor_first:  0.91  # Scale-factor of the first snaphot (cosmological run)
   time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.10  # Time difference between consecutive outputs (in internal units)
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
diff --git a/examples/EAGLE_25/run.sh b/examples/EAGLE_25/run.sh
index ae36e0d8e6ae6377682f259f506f720fbd585b29..0b6cf77d7b2461864fc24055811ee00c7dd00613 100755
--- a/examples/EAGLE_25/run.sh
+++ b/examples/EAGLE_25/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -t 16 eagle_25.yml 2>&1 | tee output.log
+../swift -c -s -G -S -t 16 eagle_25.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index 4f4b2d248ebc0c14a8b88db34de75c6f038e2f58..04c157fa86fc25f90a952e0c216285aa2235cb72 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
@@ -28,9 +28,9 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  scale_factor_first:  0.92  # Scale-factor of the first snaphot (cosmological run)
+  scale_factor_first:  0.91  # Scale-factor of the first snaphot (cosmological run)
   time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.10  # Time difference between consecutive outputs (in internal units)
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
diff --git a/examples/EAGLE_50/run.sh b/examples/EAGLE_50/run.sh
index c3b93e84192ff58e445d4e164a4d6f66d19f85fe..a0d5dee11dc58e8d19d4d0e551c5ad8eceb90548 100755
--- a/examples/EAGLE_50/run.sh
+++ b/examples/EAGLE_50/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -t 16 eagle_50.yml 2>&1 | tee output.log
+../swift -c -s -G -S -t 16 eagle_50.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index f9329eda7d8e911983018abab7fc2bda0449259e..03bb82b48394525c12b7a8df0e5c5df245fa2669 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
@@ -28,9 +28,9 @@ TimeIntegration:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  scale_factor_first:  0.92  # Scale-factor of the first snaphot (cosmological run)
+  scale_factor_first:  0.91  # Scale-factor of the first snaphot (cosmological run)
   time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.10  # Time difference between consecutive outputs (in internal units)
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
   compression:         1
 
 # Parameters governing the conserved quantities statistics
diff --git a/examples/EAGLE_6/run.sh b/examples/EAGLE_6/run.sh
index 3e1d21954d8402c9e1929119db50708859332b6c..7ef3fc2abdd1bb3fed1a228bf993bf09fb13f42c 100755
--- a/examples/EAGLE_6/run.sh
+++ b/examples/EAGLE_6/run.sh
@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -t 16 eagle_6.yml 2>&1 | tee output.log
+../swift -c -s -G -S -t 16 eagle_6.yml 2>&1 | tee output.log
 
diff --git a/examples/EAGLE_DMO_100/eagle_100.yml b/examples/EAGLE_DMO_100/eagle_100.yml
index c1bdbc2c3837ffc974c04b46ddcee9fb7bdc18c9..f04c32c8d08b5548c2c710cf8782b39a59c3821e 100644
--- a/examples/EAGLE_DMO_100/eagle_100.yml
+++ b/examples/EAGLE_DMO_100/eagle_100.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
diff --git a/examples/EAGLE_DMO_12/eagle_12.yml b/examples/EAGLE_DMO_12/eagle_12.yml
index ae4cfc64f3f4e8b046ea7624cc1e30caca23658d..2354216a5b0dcefe139d6e39699b4c67035a4173 100644
--- a/examples/EAGLE_DMO_12/eagle_12.yml
+++ b/examples/EAGLE_DMO_12/eagle_12.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
diff --git a/examples/EAGLE_DMO_25/eagle_25.yml b/examples/EAGLE_DMO_25/eagle_25.yml
index 9036827a028abd25f8afa25219e4186012dec76e..b02f9742a597687d2742b7c2d9eddf836258b06a 100644
--- a/examples/EAGLE_DMO_25/eagle_25.yml
+++ b/examples/EAGLE_DMO_25/eagle_25.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
diff --git a/examples/EAGLE_DMO_50/eagle_50.yml b/examples/EAGLE_DMO_50/eagle_50.yml
index 27bab175db4ad1e9012c9543dab0e66e41c0efbf..97299df063cd1f611f59a56ccd9b091b1217bef3 100644
--- a/examples/EAGLE_DMO_50/eagle_50.yml
+++ b/examples/EAGLE_DMO_50/eagle_50.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index 79dcd6458489958622ed7d6167d7f5f214ebd8b9..de05a9ff3c10afa7871ebeafbf4d8d272056d39f 100644
--- a/examples/ExternalPointMass/externalPointMass.yml
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
   UnitMass_in_cgs:     1.98848e33    # M_sun
-  UnitLength_in_cgs:   3.0856776e21  # kpc
+  UnitLength_in_cgs:   3.08567758e21 # kpc
   UnitVelocity_in_cgs: 1e5           # km/s
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
diff --git a/examples/ExternalPointMass/makeIC.py b/examples/ExternalPointMass/makeIC.py
index f8310e1c04265b80c6dbf6e2b003e78d95d75281..fdc5b1fd67ffcbd85beae3a9d6d1274d3d48c279 100644
--- a/examples/ExternalPointMass/makeIC.py
+++ b/examples/ExternalPointMass/makeIC.py
@@ -29,7 +29,7 @@ import random
 # physical constants in cgs
 NEWTON_GRAVITY_CGS  = 6.67408e-8
 SOLAR_MASS_IN_CGS   = 1.98848e33
-PARSEC_IN_CGS       = 3.0856776e18
+PARSEC_IN_CGS       = 3.08567758e18
 
 # choice of units
 const_unit_length_in_cgs   =   (1000*PARSEC_IN_CGS)
diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index c37a5314976c3116eb51040b1feeaa9b23ac1326..0cc11d0d8708b518b8b0b3a8df1374b6a5ead7e2 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
   UnitMass_in_cgs:     1.98848e39    # 10^6 solar masses
-  UnitLength_in_cgs:   3.0856776e21  # Kiloparsecs
+  UnitLength_in_cgs:   3.08567758e21 # Kiloparsecs
   UnitVelocity_in_cgs: 1e5           # Kilometres per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
diff --git a/examples/IsothermalPotential/isothermal.yml b/examples/IsothermalPotential/isothermal.yml
index 5dd0b831c839b3307a1c118d9bd64bda29da487c..5f626ff72e979ad0f3d404e01002be6b6018c758 100644
--- a/examples/IsothermalPotential/isothermal.yml
+++ b/examples/IsothermalPotential/isothermal.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
   UnitMass_in_cgs:     1.98848e33    # M_sun
-  UnitLength_in_cgs:   3.0856776e21  # kpc
+  UnitLength_in_cgs:   3.08567758e21 # kpc
   UnitVelocity_in_cgs: 1e5           # km/s
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
diff --git a/examples/IsothermalPotential/makeIC.py b/examples/IsothermalPotential/makeIC.py
index 85437442c6271589e3105c9c99477377b75b8d4b..eab16d21e6a4abd077dc0f4a015a4577427a3591 100644
--- a/examples/IsothermalPotential/makeIC.py
+++ b/examples/IsothermalPotential/makeIC.py
@@ -32,9 +32,8 @@ import random
 # physical constants in cgs
 NEWTON_GRAVITY_CGS  = 6.67408e-8
 SOLAR_MASS_IN_CGS   = 1.98848e33
-PARSEC_IN_CGS       = 3.0856776e18
-PROTON_MASS_IN_CGS  = 1.672621898e24
-YEAR_IN_CGS         = 3.154e+7
+PARSEC_IN_CGS       = 3.08567758e18
+YEAR_IN_CGS         = 3.15569252e7
 
 # choice of units
 const_unit_length_in_cgs   =   (1000*PARSEC_IN_CGS)
diff --git a/examples/KeplerianRing/make_movie.py b/examples/KeplerianRing/make_movie.py
index 83e783924086377a0b2aa384d2c4634bfd40443f..2e72cf143853631800b93d13abfe8633e1197380 100644
--- a/examples/KeplerianRing/make_movie.py
+++ b/examples/KeplerianRing/make_movie.py
@@ -109,7 +109,7 @@ def rotation_velocity_at_r(r, params):
         )
 
     central_mass = float(params["PointMassPotential:mass"])
-    G = 6.674e-8
+    G = 6.67408e-8
 
     v = np.sqrt( G * central_mass / r)
 
diff --git a/examples/SmallCosmoVolume/small_cosmo_volume.yml b/examples/SmallCosmoVolume/small_cosmo_volume.yml
index 3d32e0f062450f52161950a80dbf001c73b3c11d..b0057d536677f502c5c1e84a5af923042018a33b 100644
--- a/examples/SmallCosmoVolume/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/small_cosmo_volume.yml
@@ -1,7 +1,7 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
index 03445aa666c4a13f400b2458f6a2fff471248ace..481432d5875470aa464f69d5aa47fb76328cde7d 100644
--- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -1,11 +1,19 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43  # 10^10 M_su in grams
-  UnitLength_in_cgs:   3.085678e24  # Mpc in centimeters
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5   # km/s in centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
+Cosmology:
+  Omega_m: 1.
+  Omega_lambda: 0.
+  Omega_b: 1.
+  h: 1.
+  a_begin: 0.00990099
+  a_end: 1.0
+
 # Parameters governing the time integration
 TimeIntegration:
   dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
@@ -20,11 +28,11 @@ Snapshots:
 
 # Parameters governing the conserved quantities statistics
 Statistics:
-  delta_time:          10. # Time between statistics output
+  delta_time:          1.02 # Time between statistics output
 
 # Parameters for the hydrodynamics scheme
 SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation 
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
 
 # Parameters related to the initial conditions
@@ -36,14 +44,6 @@ Scheduler:
   cell_split_size:     50
   tasks_per_cell:      125
   
-Cosmology:
-  Omega_m: 1.
-  Omega_lambda: 0.
-  Omega_b: 1.
-  h: 1.
-  a_begin: 0.00990099
-  a_end: 1.0
-
 Gravity:
   mesh_side_length:   16
   eta: 0.025
diff --git a/examples/main.c b/examples/main.c
index 234bd725473486da732dd14a18fee37795357202..ba36bd87032eaa04143c1f2e147cc11a6e96875c 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1093,8 +1093,8 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
 
     fprintf(e.file_timesteps,
-            "  %6d %14e %14e %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",
-            e.step, e.time, e.cosmology->a, e.time_step, e.min_active_bin,
+            "  %6d %14e %14e %10.5f %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",
+            e.step, e.time, e.cosmology->a, e.cosmology->z, e.time_step, e.min_active_bin,
             e.max_active_bin, e.updates, e.g_updates, e.s_updates,
             e.wallclock_time, e.step_props);
     fflush(e.file_timesteps);
diff --git a/src/Makefile.am b/src/Makefile.am
index 9b043c829a20e99306eac7fa7b352859248ab7ae..fe17d77a874f3b2780b1ae00ba7523c3e63f51e6 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -64,12 +64,14 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
-		 kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h  \
+		 gravity_iact.h kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h  \
                  runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
 		 dimension.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
 		 gravity.h gravity_io.h gravity_cache.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
+		 gravity/Potential/gravity.h gravity/Potential/gravity_iact.h gravity/Potential/gravity_io.h \
+		 gravity/Potential/gravity_debug.h gravity/Potential/gravity_part.h  \
 		 sourceterms.h \
 		 equation_of_state.h \
 		 equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
diff --git a/src/cell.c b/src/cell.c
index 079950c136a569fb50d0887e079d70222ee51b00..85f8531c261ed8878bff4e32ba2419616b754372 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -264,7 +264,6 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
       temp->depth = c->depth + 1;
       temp->split = 0;
       temp->dx_max_part = 0.f;
-      temp->dx_max_gpart = 0.f;
       temp->dx_max_sort = 0.f;
       temp->nodeID = c->nodeID;
       temp->parent = c;
@@ -302,7 +301,6 @@ int cell_pack_end_step(struct cell *restrict c,
   pcells[0].ti_gravity_end_min = c->ti_gravity_end_min;
   pcells[0].ti_gravity_end_max = c->ti_gravity_end_max;
   pcells[0].dx_max_part = c->dx_max_part;
-  pcells[0].dx_max_gpart = c->dx_max_gpart;
 
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
@@ -339,7 +337,6 @@ int cell_unpack_end_step(struct cell *restrict c,
   c->ti_gravity_end_min = pcells[0].ti_gravity_end_min;
   c->ti_gravity_end_max = pcells[0].ti_gravity_end_max;
   c->dx_max_part = pcells[0].dx_max_part;
-  c->dx_max_gpart = pcells[0].dx_max_gpart;
 
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
@@ -1456,7 +1453,7 @@ void cell_activate_drift_gpart(struct cell *c, struct scheduler *s) {
         if (parent->drift_gpart == NULL)
           error("Trying to activate un-existing parent->drift_gpart");
 #endif
-	scheduler_activate(s, parent->drift_gpart);
+        scheduler_activate(s, parent->drift_gpart);
         break;
       }
     }
@@ -2609,8 +2606,6 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
   struct gpart *const gparts = c->gparts;
   struct spart *const sparts = c->sparts;
 
-  float dx_max = 0.f, dx2_max = 0.f;
-
   /* Drift irrespective of cell flags? */
   force |= c->do_grav_drift;
 
@@ -2632,15 +2627,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
 
         /* Recurse */
         cell_drift_gpart(cp, e, force);
-
-        /* Update */
-        dx_max = max(dx_max, cp->dx_max_gpart);
       }
     }
 
-    /* Store the values */
-    c->dx_max_gpart = dx_max;
-
     /* Update the time of the last drift */
     c->ti_old_gpart = ti_current;
 
@@ -2664,12 +2653,6 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
       /* Drift... */
       drift_gpart(gp, dt_drift, ti_old_gpart, ti_current);
 
-      /* Compute (square of) motion since last cell construction */
-      const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
-                        gp->x_diff[1] * gp->x_diff[1] +
-                        gp->x_diff[2] * gp->x_diff[2];
-      dx2_max = max(dx2_max, dx2);
-
       /* Init gravity force fields. */
       if (gpart_is_active(gp, e)) {
         gravity_init_gpart(gp);
@@ -2689,12 +2672,6 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
       /* Note: no need to compute dx_max as all spart have a gpart */
     }
 
-    /* Now, get the maximal particle motion from its square */
-    dx_max = sqrtf(dx2_max);
-
-    /* Store the values */
-    c->dx_max_gpart = dx_max;
-
     /* Update the time of the last drift */
     c->ti_old_gpart = ti_current;
   }
@@ -2729,8 +2706,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
     dt_drift = (ti_current - ti_old_multipole) * e->time_base;
 
   /* Drift the multipole */
-  if (ti_current > ti_old_multipole)
-    gravity_drift(c->multipole, dt_drift, c->dx_max_gpart);
+  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt_drift);
 
   /* Are we not in a leaf ? */
   if (c->split) {
@@ -2771,8 +2747,7 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
   else
     dt_drift = (ti_current - ti_old_multipole) * e->time_base;
 
-  if (ti_current > ti_old_multipole)
-    gravity_drift(c->multipole, dt_drift, c->dx_max_gpart);
+  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt_drift);
 
   /* Update the time of the last drift */
   c->ti_old_multipole = ti_current;
diff --git a/src/cell.h b/src/cell.h
index 2fcd3e20ee3d18c1b6f012f9cd4f61f364348b9f..31d525b02b49463563b21bf5aa904a8b4301f989 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -148,9 +148,6 @@ struct pcell_step {
 
   /*! Maximal distance any #part has travelled since last rebuild */
   float dx_max_part;
-
-  /*! Maximal distance any #gpart has travelled since last rebuild */
-  float dx_max_gpart;
 };
 
 /**
@@ -359,9 +356,6 @@ struct cell {
   /*! Maximum part movement in this cell since last construction. */
   float dx_max_part;
 
-  /*! Maximum gpart movement in this cell since last construction. */
-  float dx_max_gpart;
-
   /*! Nr of #part in this cell. */
   int count;
 
diff --git a/src/debug.c b/src/debug.c
index 05c21de0a73bba3a5e867a4265de0a5c14736a14..da8ef0e118b57a6aa94577898b03bcf7c56b006a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -62,7 +62,14 @@
 #error "Invalid choice of SPH variant"
 #endif
 
+/* Import the right gravity definition */
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity_debug.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_debug.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
 
 /**
  * @brief Looks for the particle with the given id and prints its information to
diff --git a/src/drift.h b/src/drift.h
index b7d5e3abe648e327f1560641676685b2d038ce3b..ff0fea744012b7143afed2a05b286d4646cdd69a 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -56,11 +56,6 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
   gp->x[0] += gp->v_full[0] * dt_drift;
   gp->x[1] += gp->v_full[1] * dt_drift;
   gp->x[2] += gp->v_full[2] * dt_drift;
-
-  /* Compute offset since last cell construction */
-  gp->x_diff[0] -= gp->v_full[0] * dt_drift;
-  gp->x_diff[1] -= gp->v_full[1] * dt_drift;
-  gp->x_diff[2] -= gp->v_full[2] * dt_drift;
 }
 
 /**
diff --git a/src/engine.c b/src/engine.c
index 4ad4d31d87b31907928e0c35a0c2a34f2cc29584..52b252e2c16f1f5f392db26fa69bc63d1ba7afe9 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -346,7 +346,7 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
                                            task_subtype_none, 0, 0, c, NULL);
 
         if (periodic) scheduler_addunlock(s, c->drift_gpart, c->grav_mesh);
-        if (periodic) scheduler_addunlock(s, c->grav_mesh, c->super->end_force);
+        if (periodic) scheduler_addunlock(s, c->grav_mesh, c->grav_down);
         scheduler_addunlock(s, c->init_grav, c->grav_long_range);
         scheduler_addunlock(s, c->grav_long_range, c->grav_down);
         scheduler_addunlock(s, c->grav_down, c->super->end_force);
@@ -2378,8 +2378,8 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
   int delta_p = delta;
 
   /* Special case where every cell is in range of every other one */
-  if(delta >= cdim[0] / 2) {
-    if(cdim[0] % 2 == 0) {
+  if (delta >= cdim[0] / 2) {
+    if (cdim[0] % 2 == 0) {
       delta_m = cdim[0] / 2;
       delta_p = cdim[0] / 2 - 1;
     } else {
@@ -2662,13 +2662,14 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
     /* Link self tasks to cells. */
     else if (t_type == task_type_self) {
       atomic_inc(&ci->nr_tasks);
+
       if (t_subtype == task_subtype_density) {
         engine_addlink(e, &ci->density, t);
       }
-      if (t_subtype == task_subtype_grav) {
+      else if (t_subtype == task_subtype_grav) {
         engine_addlink(e, &ci->grav, t);
       }
-      if (t_subtype == task_subtype_external_grav) {
+      else if (t_subtype == task_subtype_external_grav) {
         engine_addlink(e, &ci->grav, t);
       }
 
@@ -2676,16 +2677,17 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
     } else if (t_type == task_type_pair) {
       atomic_inc(&ci->nr_tasks);
       atomic_inc(&cj->nr_tasks);
+
       if (t_subtype == task_subtype_density) {
         engine_addlink(e, &ci->density, t);
         engine_addlink(e, &cj->density, t);
       }
-      if (t_subtype == task_subtype_grav) {
+      else if (t_subtype == task_subtype_grav) {
         engine_addlink(e, &ci->grav, t);
         engine_addlink(e, &cj->grav, t);
       }
 #ifdef SWIFT_DEBUG_CHECKS
-      if (t_subtype == task_subtype_external_grav) {
+      else if (t_subtype == task_subtype_external_grav) {
         error("Found a pair/external-gravity task...");
       }
 #endif
@@ -2693,13 +2695,14 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
       /* Link sub-self tasks to cells. */
     } else if (t_type == task_type_sub_self) {
       atomic_inc(&ci->nr_tasks);
+
       if (t_subtype == task_subtype_density) {
         engine_addlink(e, &ci->density, t);
       }
-      if (t_subtype == task_subtype_grav) {
+      else if (t_subtype == task_subtype_grav) {
         engine_addlink(e, &ci->grav, t);
       }
-      if (t_subtype == task_subtype_external_grav) {
+      else if (t_subtype == task_subtype_external_grav) {
         engine_addlink(e, &ci->grav, t);
       }
 
@@ -2707,16 +2710,17 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
     } else if (t_type == task_type_sub_pair) {
       atomic_inc(&ci->nr_tasks);
       atomic_inc(&cj->nr_tasks);
+
       if (t_subtype == task_subtype_density) {
         engine_addlink(e, &ci->density, t);
         engine_addlink(e, &cj->density, t);
       }
-      if (t_subtype == task_subtype_grav) {
+      else if (t_subtype == task_subtype_grav) {
         engine_addlink(e, &ci->grav, t);
         engine_addlink(e, &cj->grav, t);
       }
 #ifdef SWIFT_DEBUG_CHECKS
-      if (t_subtype == task_subtype_external_grav) {
+      else if (t_subtype == task_subtype_external_grav) {
         error("Found a sub-pair/external-gravity task...");
       }
 #endif
@@ -3188,15 +3192,26 @@ void engine_maketasks(struct engine *e) {
   /* Re-set the scheduler. */
   scheduler_reset(sched, engine_estimate_nr_tasks(e));
 
+  ticks tic2 = getticks();
+
   /* Construct the firt hydro loop over neighbours */
-  if (e->policy & engine_policy_hydro) {
+  if (e->policy & engine_policy_hydro)
     threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper, NULL,
                    s->nr_cells, 1, 0, e);
-  }
+
+  if (e->verbose)
+    message("Making hydro tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
+  tic2 = getticks();
 
   /* Add the self gravity tasks. */
   if (e->policy & engine_policy_self_gravity) engine_make_self_gravity_tasks(e);
 
+  if (e->verbose)
+    message("Making gravity tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
   /* Add the external gravity tasks. */
   if (e->policy & engine_policy_external_gravity)
     engine_make_external_gravity_tasks(e);
@@ -3234,9 +3249,15 @@ void engine_maketasks(struct engine *e) {
     error("Failed to allocate cell-task links.");
   e->nr_links = 0;
 
+  tic2 = getticks();
+
   /* Split the tasks. */
   scheduler_splittasks(sched);
 
+  if (e->verbose)
+    message("Splitting tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that we are not left with invalid tasks */
   for (int i = 0; i < e->sched.nr_tasks; ++i) {
@@ -3245,21 +3266,35 @@ void engine_maketasks(struct engine *e) {
   }
 #endif
 
+  tic2 = getticks();
+
   /* Count the number of tasks associated with each cell and
      store the density tasks in each cell, and make each sort
      depend on the sorts of its progeny. */
   threadpool_map(&e->threadpool, engine_count_and_link_tasks_mapper,
                  sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
 
+  if (e->verbose)
+    message("Counting and linking tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
+  tic2 = getticks();
+
   /* Now that the self/pair tasks are at the right level, set the super
    * pointers. */
   threadpool_map(&e->threadpool, cell_set_super_mapper, cells, nr_cells,
                  sizeof(struct cell), 0, e);
 
+  if (e->verbose)
+    message("Setting super-pointers took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
   /* Append hierarchical tasks to each cell. */
   threadpool_map(&e->threadpool, engine_make_hierarchical_tasks_mapper, cells,
                  nr_cells, sizeof(struct cell), 0, e);
 
+  tic2 = getticks();
+
   /* Run through the tasks and make force tasks for each density task.
      Each force task depends on the cell ghosts and unlocks the kick task
      of its super-cell. */
@@ -3267,10 +3302,20 @@ void engine_maketasks(struct engine *e) {
     threadpool_map(&e->threadpool, engine_make_extra_hydroloop_tasks_mapper,
                    sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
 
+  if (e->verbose)
+    message("Making extra hydroloop tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
+  tic2 = getticks();
+
   /* Add the dependencies for the gravity stuff */
   if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
     engine_link_gravity_tasks(e);
 
+  if (e->verbose)
+    message("Linking gravity tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
 #ifdef WITH_MPI
 
   /* Add the communication tasks if MPI is being used. */
@@ -3323,15 +3368,33 @@ void engine_maketasks(struct engine *e) {
   }
 #endif
 
+  tic2 = getticks();
+
   /* Set the unlocks per task. */
   scheduler_set_unlocks(sched);
 
+  if (e->verbose)
+    message("Setting unlocks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
+  tic2 = getticks();
+
   /* Rank the tasks. */
   scheduler_ranktasks(sched);
 
+  if (e->verbose)
+    message("Ranking the tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
+  tic2 = getticks();
+
   /* Weight the tasks. */
   scheduler_reweight(sched, e->verbose);
 
+  if (e->verbose)
+    message("Reweighting tasks took %.3f %s (including reweight).",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
   /* Set the tasks age. */
   e->tasks_age = 0;
 
@@ -3899,13 +3962,14 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
 
   /* Clear the forcerebuild flag, whatever it was. */
   e->forcerebuild = 0;
+  e->restarting = 0;
 
   /* Re-build the space. */
   space_rebuild(e->s, e->verbose);
 
   /* Re-compute the mesh forces */
   if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
-    pm_mesh_compute_potential(e->mesh, e);
+    pm_mesh_compute_potential(e->mesh, e->s, e->verbose);
 
   /* Re-compute the maximal RMS displacement constraint */
   if (e->policy & engine_policy_cosmology)
@@ -3975,7 +4039,7 @@ void engine_prepare(struct engine *e) {
   const ticks tic = getticks();
 
   /* Unskip active tasks and check for rebuild */
-  if (!e->forcerepart) engine_unskip(e);
+  if (!e->forcerepart && !e->restarting) engine_unskip(e);
 
 #ifdef WITH_MPI
   MPI_Allreduce(MPI_IN_PLACE, &e->forcerebuild, 1, MPI_INT, MPI_MAX,
@@ -3989,7 +4053,7 @@ void engine_prepare(struct engine *e) {
   if (e->forcerebuild) {
 
     /* Let's start by drifting everybody to the current time */
-    engine_drift_all(e);
+    if (!e->restarting) engine_drift_all(e);
 
     /* And rebuild */
     engine_rebuild(e, 0);
@@ -4447,6 +4511,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   struct clocks_time time1, time2;
   clocks_gettime(&time1);
 
+  /* Update the softening lengths */
+  if (e->policy & engine_policy_self_gravity)
+    gravity_update(e->gravity_properties, e->cosmology);
+
   /* Start by setting the particles in a good state */
   if (e->nodeID == 0) message("Setting particles to a valid state...");
   engine_first_init_particles(e);
@@ -4640,20 +4708,25 @@ void engine_step(struct engine *e) {
 
     /* Print some information to the screen */
     printf(
-        "  %6d %14e %14e %10.5f %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",
+	"  %6d %14e %14e %10.5f %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",
         e->step, e->time, e->cosmology->a, e->cosmology->z, e->time_step,
         e->min_active_bin, e->max_active_bin, e->updates, e->g_updates,
         e->s_updates, e->wallclock_time, e->step_props);
     fflush(stdout);
 
-    fprintf(e->file_timesteps,
-            "  %6d %14e %14e %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",
-            e->step, e->time, e->cosmology->a, e->time_step, e->min_active_bin,
-            e->max_active_bin, e->updates, e->g_updates, e->s_updates,
-            e->wallclock_time, e->step_props);
+    if(!e->restarting)
+      fprintf(e->file_timesteps,
+	      "  %6d %14e %14e %10.5f %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",
+	      e->step, e->time, e->cosmology->a, e->cosmology->z, e->time_step, e->min_active_bin,
+	      e->max_active_bin, e->updates, e->g_updates, e->s_updates,
+	      e->wallclock_time, e->step_props);
     fflush(e->file_timesteps);
   }
 
+  /* We need some cells to exist but not the whole task stuff. */
+  if(e->restarting)
+    space_rebuild(e->s, e->verbose);
+
   /* Move forward in time */
   e->ti_old = e->ti_current;
   e->ti_current = e->ti_end_min;
@@ -4662,6 +4735,11 @@ void engine_step(struct engine *e) {
   e->step += 1;
   e->step_props = engine_step_prop_none;
 
+  /* When restarting, move everyone to the current time. */
+  if(e->restarting)
+    engine_drift_all(e);
+
+  /* Get the physical value of the time and time-step size */
   if (e->policy & engine_policy_cosmology) {
     e->time_old = e->time;
     cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
@@ -4923,6 +5001,10 @@ void engine_dump_restarts(struct engine *e, int drifted_all, int force) {
     MPI_Bcast(&dump, 1, MPI_INT, 0, MPI_COMM_WORLD);
 #endif
     if (dump) {
+
+      if(e->nodeID == 0)
+	message("Writing restart files");
+
       /* Clean out the previous saved files, if found. Do this now as we are
        * MPI synchronized. */
       restart_remove_previous(e->restart_file);
@@ -5758,6 +5840,7 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
   e->nr_proxies = 0;
   e->forcerebuild = 1;
   e->forcerepart = 0;
+  e->restarting = restart;
   e->step_props = engine_step_prop_none;
   e->links = NULL;
   e->nr_links = 0;
@@ -5968,8 +6051,8 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
               engine_step_prop_snapshot, engine_step_prop_restarts);
 
       fprintf(e->file_timesteps,
-              "# %6s %14s %14s %14s %9s %12s %12s %12s %16s [%s] %6s\n", "Step",
-              "Time", "Scale-factor", "Time-step", "Time-bins", "Updates",
+              "# %6s %14s %14s %10s %14s %9s %12s %12s %12s %16s [%s] %6s\n", "Step",
+              "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins", "Updates",
               "g-Updates", "s-Updates", "Wall-clock time", clocks_getunit(),
               "Props");
       fflush(e->file_timesteps);
@@ -6022,9 +6105,10 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
         "dt=%e",
         e->time_base);
 
-  if (e->dt_max > (e->time_end - e->time_begin) && e->nodeID == 0)
-    error("Maximal time-step size larger than the simulation run time t=%e",
-          e->time_end - e->time_begin);
+  if (!(e->policy & engine_policy_cosmology))
+    if (e->dt_max > (e->time_end - e->time_begin) && e->nodeID == 0)
+      error("Maximal time-step size larger than the simulation run time t=%e",
+            e->time_end - e->time_begin);
 
   /* Deal with outputs */
   if (e->policy & engine_policy_cosmology) {
@@ -6070,6 +6154,17 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
           e->time_first_statistics, e->time_begin);
   }
 
+  /* Get the total mass */
+  e->total_mass = 0.;
+  for (size_t i = 0; i < e->s->nr_gparts; ++i)
+    e->total_mass += e->s->gparts[i].mass;
+
+/* Reduce the total mass */
+#ifdef WITH_MPI
+  MPI_Allreduce(MPI_IN_PLACE, &e->total_mass, 1, MPI_DOUBLE, MPI_SUM,
+                MPI_COMM_WORLD);
+#endif
+
   /* Find the time of the first snapshot  output */
   engine_compute_next_snapshot_time(e);
 
diff --git a/src/engine.h b/src/engine.h
index b29b5bb0bb6013e14c49f0d45bd072984822e012..ae430d9f796d6cf9bd4a1973f15b3a4e986b9997 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -200,6 +200,9 @@ struct engine {
   /* Total numbers of particles in the system. */
   long long total_nr_parts, total_nr_gparts, total_nr_sparts;
 
+  /* Total mass in the simulation */
+  double total_mass;
+
   /* The internal system of units */
   const struct unit_system *internal_units;
 
@@ -261,6 +264,9 @@ struct engine {
   /* Wallclock time of the last time-step */
   float wallclock_time;
 
+  /* Are we in the process of restaring a simulation? */
+  int restarting;
+
   /* Force the engine to rebuild? */
   int forcerebuild;
 
diff --git a/src/gravity.c b/src/gravity.c
index d0eb589d15ba89aecc5a581de0c143dc7504d3ee..24aa7f56b9a950227b04a92dde1ec5d9af4a84ed 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -152,6 +152,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
     bzero(fewald_z, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
     bzero(potewald, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
 
+    /* Hernquist, Bouchet & Suto, 1991, Eq. 2.10 and just below Eq. 2.15 */
     potewald[0][0][0] = 2.8372975f;
 
     /* Compute the values in one of the octants */
diff --git a/src/gravity.h b/src/gravity.h
index 6497de8294dfa3f207332ff696ddb992c875eb28..57cad8eba5f0772f85affd031a450c3ecb4dde59 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -27,10 +27,16 @@
 #include "inline.h"
 #include "part.h"
 
-/* So far only one model here */
-/* Straight-forward import */
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity.h"
-#include "./gravity/Default/gravity_iact.h"
+#define GRAVITY_IMPLEMENTATION "Default (no potential)"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity.h"
+#define GRAVITY_IMPLEMENTATION "With potential calculation"
+#else
+#error "Invalid choice of gravity variant"
+#endif
 
 struct engine;
 struct space;
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 9f0db3f3bde9aef7a847d22dbc36b35b192b9304..2713c9ee7affca4f06b369d038916f76b8c2ee48 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -51,19 +51,36 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
 }
 
 /**
- * @brief Returns the comoving potential of a particle
+ * @brief Add a contribution to this particle's potential.
+ *
+ * Here we do nothing as this version does not accumulate potential.
+ *
+ * @param gp The particle.
+ * @param pot The contribution to add.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {}
+
+/**
+ * @brief Returns the comoving potential of a particle.
+ *
+ * This returns 0 as this flavour of gravity does not compute the
+ * particles' potential.
  *
  * @param gp The particle of interest
  */
 __attribute__((always_inline)) INLINE static float
 gravity_get_comoving_potential(const struct gpart* restrict gp) {
 
-  return gp->potential;
+  return 0.f;
 }
 
 /**
  * @brief Returns the physical potential of a particle
  *
+ * This returns 0 as this flavour of gravity does not compute the
+ * particles' potential.
+ *
  * @param gp The particle of interest.
  * @param cosmo The cosmological model.
  */
@@ -71,7 +88,7 @@ __attribute__((always_inline)) INLINE static float
 gravity_get_physical_potential(const struct gpart* restrict gp,
                                const struct cosmology* cosmo) {
 
-  return gp->potential * cosmo->a_inv;
+  return 0.f;
 }
 
 /**
@@ -128,7 +145,6 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
   gp->a_grav[0] = 0.f;
   gp->a_grav[1] = 0.f;
   gp->a_grav[2] = 0.f;
-  gp->potential = 0.f;
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = 0.f;
@@ -145,19 +161,25 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
 /**
  * @brief Finishes the gravity calculation.
  *
- * Multiplies the forces and accelerations by the appropiate constants
+ * Multiplies the forces and accelerations by the appropiate constants.
+ * Applies cosmological correction for periodic BCs.
+ *
+ * No need to apply the potential normalisation correction for periodic
+ * BCs here since we do not compute the potential.
  *
  * @param gp The particle to act upon
- * @param const_G Newton's constant in internal units
+ * @param const_G Newton's constant in internal units.
+ * @param potential_normalisation Term to be added to all the particles.
+ * @param periodic Are we using periodic BCs?
  */
 __attribute__((always_inline)) INLINE static void gravity_end_force(
-    struct gpart* gp, float const_G) {
+    struct gpart* gp, float const_G, const float potential_normalisation,
+    const int periodic) {
 
   /* Let's get physical... */
   gp->a_grav[0] *= const_G;
   gp->a_grav[1] *= const_G;
   gp->a_grav[2] *= const_G;
-  gp->potential *= const_G;
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM *= const_G;
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index ad477c54f49f9d30ee492d84012d7a9228401c5f..6fd7ae205ecb124344f0780a5e8a276445ad07b4 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -55,21 +55,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
 
     /* Get Newtonian gravity */
     *f_ij = mass * r_inv * r_inv * r_inv;
-    *pot_ij = -mass * r_inv;
 
   } else {
 
     const float r = r2 * r_inv;
     const float ui = r * h_inv;
 
-    float W_f_ij, W_pot_ij;
+    float W_f_ij;
     kernel_grav_force_eval(ui, &W_f_ij);
-    kernel_grav_pot_eval(ui, &W_pot_ij);
 
     /* Get softened gravity */
     *f_ij = mass * h_inv3 * W_f_ij;
-    *pot_ij = mass * h_inv * W_pot_ij;
   }
+
+  /* No potential calculation */
+  *pot_ij = 0.f;
 }
 
 /**
@@ -101,28 +101,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
 
     /* Get Newtonian gravity */
     *f_ij = mass * r_inv * r_inv * r_inv;
-    *pot_ij = -mass * r_inv;
 
   } else {
 
     const float ui = r * h_inv;
-    float W_f_ij, W_pot_ij;
+    float W_f_ij;
 
     kernel_grav_force_eval(ui, &W_f_ij);
-    kernel_grav_pot_eval(ui, &W_pot_ij);
 
     /* Get softened gravity */
     *f_ij = mass * h_inv3 * W_f_ij;
-    *pot_ij = mass * h_inv * W_pot_ij;
   }
 
   /* Get long-range correction */
   const float u_lr = r * r_s_inv;
-  float corr_f_lr, corr_pot_lr;
+  float corr_f_lr;
   kernel_long_grav_force_eval(u_lr, &corr_f_lr);
-  kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
   *f_ij *= corr_f_lr;
-  *pot_ij *= corr_pot_lr;
+
+  /* No potential calculation */
+  *pot_ij = 0.f;
 }
 
 /**
@@ -160,7 +158,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_x = f_ij * r_x;
   *f_y = f_ij * r_y;
   *f_z = f_ij * r_z;
-  *pot = pot_ij;
 
 #else
 
@@ -176,7 +173,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_x = m->M_000 * d.D_100;
   *f_y = m->M_000 * d.D_010;
   *f_z = m->M_000 * d.D_001;
-  *pot = m->M_000 * d.D_000;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
@@ -187,7 +183,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
 /* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
 /* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
 /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
-/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
 
 #endif
 
@@ -200,8 +195,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
   *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
           m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
-  *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
-          m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
 
 #endif
 
@@ -220,10 +213,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
           m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
           m->M_300 * d.D_301;
-  *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
-          m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
-          m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
-          m->M_300 * d.D_300;
 
 #endif
 
@@ -231,8 +220,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_x *= -1.f;
   *f_y *= -1.f;
   *f_z *= -1.f;
-  *pot *= -1.f;
 #endif
+
+  /* No potential calculation */
+  *pot = 0.f;
 }
 
 /**
@@ -272,7 +263,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_x = f_ij * r_x;
   *f_y = f_ij * r_y;
   *f_z = f_ij * r_z;
-  *pot = -pot_ij;
 
 #else
 
@@ -288,7 +278,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_x = m->M_000 * d.D_100;
   *f_y = m->M_000 * d.D_010;
   *f_z = m->M_000 * d.D_001;
-  *pot = m->M_000 * d.D_000;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
@@ -299,7 +288,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
 /* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
 /* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
 /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
-/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
 
 #endif
 
@@ -312,8 +300,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
   *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
           m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
-  *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
-          m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
 
 #endif
 
@@ -332,19 +318,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
           m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
           m->M_300 * d.D_301;
-  *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
-          m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
-          m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
-          m->M_300 * d.D_300;
-
 #endif
 
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
   *f_z *= -1.f;
-  *pot *= -1.f;
 #endif
+
+  /* No potential calculation */
+  *pot = 0.f;
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 7f453179641e2ba16b30e3172ddd7853245a1d2f..1ba3899e7ecc346227c10679bb8b704937c625b2 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -104,7 +104,7 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
                                               int* num_fields) {
 
   /* Say how much we want to write */
-  *num_fields = 5;
+  *num_fields = 4;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_gpart(
@@ -115,8 +115,6 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
       io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass);
   list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                  UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
-  list[4] = io_make_output_field("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL,
-                                 gparts, potential);
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index d325baf4de0938d8df539d38ae10f8f3f3ec7d2b..bd73c56da82877415f5abc9edf41ede1c551f16f 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -29,9 +29,6 @@ struct gpart {
   /*! Particle position. */
   double x[3];
 
-  /*! Offset between current position and position at last tree rebuild. */
-  float x_diff[3];
-
   /*! Particle velocity. */
   float v_full[3];
 
@@ -41,9 +38,6 @@ struct gpart {
   /*! Particle mass. */
   float mass;
 
-  /*! Gravitational potential */
-  float potential;
-
   /*! Time-step length */
   timebin_t time_bin;
 
diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a6c0fba18856b57911d49bcee6915f5003e2e68
--- /dev/null
+++ b/src/gravity/Potential/gravity.h
@@ -0,0 +1,223 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *               2016 Tom Theuns (tom.theuns@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_POTENTIAL_GRAVITY_H
+#define SWIFT_POTENTIAL_GRAVITY_H
+
+#include <float.h>
+
+#include "cosmology.h"
+#include "gravity_properties.h"
+#include "kernel_gravity.h"
+#include "minmax.h"
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_mass(
+    const struct gpart* restrict gp) {
+
+  return gp->mass;
+}
+
+/**
+ * @brief Returns the softening of a particle
+ *
+ * @param gp The particle of interest
+ * @param grav_props The global gravity properties.
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_softening(
+    const struct gpart* gp, const struct gravity_props* restrict grav_props) {
+
+  return grav_props->epsilon_cur;
+}
+
+/**
+ * @brief Add a contribution to this particle's potential.
+ *
+ * @param gp The particle.
+ * @param pot The contribution to add.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {
+
+  gp->potential += pot;
+}
+
+/**
+ * @brief Returns the comoving potential of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_comoving_potential(const struct gpart* restrict gp) {
+
+  return gp->potential;
+}
+
+/**
+ * @brief Returns the physical potential of a particle
+ *
+ * @param gp The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_physical_potential(const struct gpart* restrict gp,
+                               const struct cosmology* cosmo) {
+
+  return gp->potential * cosmo->a_inv;
+}
+
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * We use Gadget-2's type 0 time-step criterion.
+ *
+ * @param gp Pointer to the g-particle data.
+ * @param a_hydro The accelerations coming from the hydro scheme (can be 0).
+ * @param grav_props Constants used in the gravity scheme.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_self(const struct gpart* const gp,
+                              const float a_hydro[3],
+                              const struct gravity_props* restrict grav_props,
+                              const struct cosmology* cosmo) {
+
+  /* Get physical acceleration (gravity contribution) */
+  float a_phys_x = gp->a_grav[0] * cosmo->a_factor_grav_accel;
+  float a_phys_y = gp->a_grav[1] * cosmo->a_factor_grav_accel;
+  float a_phys_z = gp->a_grav[2] * cosmo->a_factor_grav_accel;
+
+  /* Get physical acceleration (hydro contribution) */
+  a_phys_x += a_hydro[0] * cosmo->a_factor_hydro_accel;
+  a_phys_y += a_hydro[1] * cosmo->a_factor_hydro_accel;
+  a_phys_z += a_hydro[2] * cosmo->a_factor_hydro_accel;
+
+  const float ac2 =
+      a_phys_x * a_phys_x + a_phys_y * a_phys_y + a_phys_z * a_phys_z;
+
+  const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
+
+  const float epsilon = gravity_get_softening(gp, grav_props);
+
+  const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
+                         cosmo->a * grav_props->eta * epsilon * ac_inv);
+
+  return dt;
+}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void gravity_init_gpart(
+    struct gpart* gp) {
+
+  /* Zero the acceleration */
+  gp->a_grav[0] = 0.f;
+  gp->a_grav[1] = 0.f;
+  gp->a_grav[2] = 0.f;
+  gp->potential = 0.f;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->potential_PM = 0.f;
+  gp->a_grav_PM[0] = 0.f;
+  gp->a_grav_PM[1] = 0.f;
+  gp->a_grav_PM[2] = 0.f;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+  gp->num_interacted = 0;
+#endif
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants.
+ * Applies cosmological correction for periodic BCs.
+ *
+ * @param gp The particle to act upon
+ * @param const_G Newton's constant in internal units.
+ * @param potential_normalisation Term to be added to all the particles.
+ */
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+    struct gpart* gp, float const_G, const float potential_normalisation,
+    const int periodic) {
+
+  /* Apply the periodic correction to the peculiar potential */
+  if (periodic) gp->potential += potential_normalisation;
+
+  /* Let's get physical... */
+  gp->a_grav[0] *= const_G;
+  gp->a_grav[1] *= const_G;
+  gp->a_grav[2] *= const_G;
+  gp->potential *= const_G;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->potential_PM *= const_G;
+  gp->a_grav_PM[0] *= const_G;
+  gp->a_grav_PM[1] *= const_G;
+  gp->a_grav_PM[2] *= const_G;
+#endif
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+    struct gpart* gp, float dt) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param gp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_reset_predicted_values(struct gpart* gp) {}
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ * @param grav_props The global properties of the gravity calculation.
+ */
+__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
+    struct gpart* gp, const struct gravity_props* grav_props) {
+
+  gp->time_bin = 0;
+
+  gravity_init_gpart(gp);
+}
+
+#endif /* SWIFT_POTENTIAL_GRAVITY_H */
diff --git a/src/gravity/Potential/gravity_debug.h b/src/gravity/Potential/gravity_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..621f3ccb9c4621f44f902c37a2e03bd7575defe7
--- /dev/null
+++ b/src/gravity/Potential/gravity_debug.h
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_POTENTIAL_GRAVITY_DEBUG_H
+#define SWIFT_POTENTIAL_GRAVITY_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void gravity_debug_particle(
+    const struct gpart* p) {
+  printf(
+      "mass=%.3e time_bin=%d\n"
+      "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e] "
+      "pot=%.5e\n",
+      p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
+      p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2],
+      p->potential);
+#ifdef SWIFT_DEBUG_CHECKS
+  printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
+         p->ti_drift, p->ti_kick);
+#endif
+}
+
+#endif /* SWIFT_POTENTIAL_GRAVITY_DEBUG_H */
diff --git a/src/gravity/Potential/gravity_iact.h b/src/gravity/Potential/gravity_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d6029b96214549686e00789dcb8a3a3ad2750b4
--- /dev/null
+++ b/src/gravity/Potential/gravity_iact.h
@@ -0,0 +1,350 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    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_POTENTIAL_GRAVITY_IACT_H
+#define SWIFT_POTENTIAL_GRAVITY_IACT_H
+
+/* Includes. */
+#include "kernel_gravity.h"
+#include "kernel_long_gravity.h"
+#include "multipole.h"
+
+/* Standard headers */
+#include <float.h>
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
+    const float r2, const float h2, const float h_inv, const float h_inv3,
+    const float mass, float *f_ij, float *pot_ij) {
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+
+  /* Should we soften ? */
+  if (r2 >= h2) {
+
+    /* Get Newtonian gravity */
+    *f_ij = mass * r_inv * r_inv * r_inv;
+    *pot_ij = -mass * r_inv;
+
+  } else {
+
+    const float r = r2 * r_inv;
+    const float ui = r * h_inv;
+
+    float W_f_ij, W_pot_ij;
+    kernel_grav_force_eval(ui, &W_f_ij);
+    kernel_grav_pot_eval(ui, &W_pot_ij);
+
+    /* Get softened gravity */
+    *f_ij = mass * h_inv3 * W_f_ij;
+    *pot_ij = mass * h_inv * W_pot_ij;
+  }
+}
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass truncated for long-distance periodicity.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param r_s_inv Inverse of the mesh smoothing scale.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
+    const float r2, const float h2, const float h_inv, const float h_inv3,
+    const float mass, const float r_s_inv, float *f_ij, float *pot_ij) {
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+  const float r = r2 * r_inv;
+
+  /* Should we soften ? */
+  if (r2 >= h2) {
+
+    /* Get Newtonian gravity */
+    *f_ij = mass * r_inv * r_inv * r_inv;
+    *pot_ij = -mass * r_inv;
+
+  } else {
+
+    const float ui = r * h_inv;
+    float W_f_ij, W_pot_ij;
+
+    kernel_grav_force_eval(ui, &W_f_ij);
+    kernel_grav_pot_eval(ui, &W_pot_ij);
+
+    /* Get softened gravity */
+    *f_ij = mass * h_inv3 * W_f_ij;
+    *pot_ij = mass * h_inv * W_pot_ij;
+  }
+
+  /* Get long-range correction */
+  const float u_lr = r * r_s_inv;
+  float corr_f_lr, corr_pot_lr;
+  kernel_long_grav_force_eval(u_lr, &corr_f_lr);
+  kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
+  *f_ij *= corr_f_lr;
+  *pot_ij *= corr_pot_lr;
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole and trace of the octupole terms only and defaults to
+ * the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
+    const float r_x, const float r_y, const float r_z, const float r2,
+    const float h, const float h_inv, const struct multipole *m, float *f_x,
+    float *f_y, float *f_z, float *pot) {
+
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
+  float f_ij, pot_ij;
+  runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
+                           &f_ij, &pot_ij);
+  *f_x = f_ij * r_x;
+  *f_y = f_ij * r_y;
+  *f_z = f_ij * r_z;
+  *pot = pot_ij;
+
+#else
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2);
+
+  /* Compute the derivatives of the potential */
+  struct potential_derivatives_M2P d;
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
+                                    &d);
+
+  /* 0th order contributions */
+  *f_x = m->M_000 * d.D_100;
+  *f_y = m->M_000 * d.D_010;
+  *f_z = m->M_000 * d.D_001;
+  *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+/* 1st order contributions */
+
+/* 1st order contributions are all 0 since the dipole is 0 */
+
+/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+  /* 2nd order contributions */
+  *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+          m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+  *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+          m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+  *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+          m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+  *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+          m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+  /* 3rd order contributions */
+  *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+          m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+          m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+          m->M_300 * d.D_400;
+  *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+          m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+          m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+          m->M_300 * d.D_310;
+  *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+          m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+          m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+          m->M_300 * d.D_301;
+  *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+          m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+          m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+          m->M_300 * d.D_300;
+
+#endif
+
+  /* Take care of the the sign convention */
+  *f_x *= -1.f;
+  *f_y *= -1.f;
+  *f_z *= -1.f;
+  *pot *= -1.f;
+#endif
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole, truncated for
+ * long-range periodicity.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole term and trace of the octupole terms only and
+ * defaults to the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param r_s_inv The inverse of the gravity mesh-smoothing scale.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
+    const float r_x, const float r_y, const float r_z, const float r2,
+    const float h, const float h_inv, const float r_s_inv,
+    const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
+
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
+  float f_ij, pot_ij;
+  runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
+                                m->M_000, r_s_inv, &f_ij, &pot_ij);
+  *f_x = f_ij * r_x;
+  *f_y = f_ij * r_y;
+  *f_z = f_ij * r_z;
+  *pot = -pot_ij;
+
+#else
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2);
+
+  /* Compute the derivatives of the potential */
+  struct potential_derivatives_M2P d;
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
+                                    r_s_inv, &d);
+
+  /* 0th order contributions */
+  *f_x = m->M_000 * d.D_100;
+  *f_y = m->M_000 * d.D_010;
+  *f_z = m->M_000 * d.D_001;
+  *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+/* 1st order contributions */
+
+/* 1st order contributions are all 0 since the dipole is 0 */
+
+/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+  /* 2nd order contributions */
+  *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+          m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+  *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+          m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+  *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+          m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+  *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+          m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+  /* 3rd order contributions */
+  *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+          m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+          m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+          m->M_300 * d.D_400;
+  *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+          m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+          m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+          m->M_300 * d.D_310;
+  *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+          m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+          m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+          m->M_300 * d.D_301;
+  *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+          m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+          m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+          m->M_300 * d.D_300;
+
+#endif
+
+  /* Take care of the the sign convention */
+  *f_x *= -1.f;
+  *f_y *= -1.f;
+  *f_z *= -1.f;
+  *pot *= -1.f;
+#endif
+}
+
+#endif /* SWIFT_POTENTIAL_GRAVITY_IACT_H */
diff --git a/src/gravity/Potential/gravity_io.h b/src/gravity/Potential/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..6aa4cbb4786af99ac372564ed67f4ce77c08f25c
--- /dev/null
+++ b/src/gravity/Potential/gravity_io.h
@@ -0,0 +1,122 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_POTENTIAL_GRAVITY_IO_H
+#define SWIFT_POTENTIAL_GRAVITY_IO_H
+
+#include "io_properties.h"
+
+INLINE static void convert_gpart_pos(const struct engine* e,
+                                     const struct gpart* gp, double* ret) {
+
+  if (e->s->periodic) {
+    ret[0] = box_wrap(gp->x[0], 0.0, e->s->dim[0]);
+    ret[1] = box_wrap(gp->x[1], 0.0, e->s->dim[1]);
+    ret[2] = box_wrap(gp->x[2], 0.0, e->s->dim[2]);
+  } else {
+    ret[0] = gp->x[0];
+    ret[1] = gp->x[1];
+    ret[2] = gp->x[2];
+  }
+}
+
+INLINE static void convert_gpart_vel(const struct engine* e,
+                                     const struct gpart* gp, float* ret) {
+
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const struct cosmology* cosmo = e->cosmology;
+  const integertime_t ti_current = e->ti_current;
+  const double time_base = e->time_base;
+
+  const integertime_t ti_beg = get_integer_time_begin(ti_current, gp->time_bin);
+  const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+
+  /* Get time-step since the last kick */
+  float dt_kick_grav;
+  if (with_cosmology) {
+    dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_grav -=
+        cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+  } else {
+    dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+  }
+
+  /* Extrapolate the velocites to the current time */
+  ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+  ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+  ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+  /* Conversion from internal units to peculiar velocities */
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
+}
+
+/**
+ * @brief Specifies which g-particle fields to read from a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void darkmatter_read_particles(struct gpart* gparts,
+                                             struct io_props* list,
+                                             int* num_fields) {
+
+  /* Say how much we want to read */
+  *num_fields = 4;
+
+  /* List what we want to read */
+  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+                                UNIT_CONV_LENGTH, gparts, x);
+  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+                                UNIT_CONV_SPEED, gparts, v_full);
+  list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+                                gparts, mass);
+  list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+                                UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
+}
+
+/**
+ * @brief Specifies which g-particle fields to write to a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void darkmatter_write_particles(const struct gpart* gparts,
+                                              struct io_props* list,
+                                              int* num_fields) {
+
+  /* Say how much we want to write */
+  *num_fields = 5;
+
+  /* List what we want to write */
+  list[0] = io_make_output_field_convert_gpart(
+      "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos);
+  list[1] = io_make_output_field_convert_gpart(
+      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, gparts, convert_gpart_vel);
+  list[2] =
+      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass);
+  list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+                                 UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
+  list[4] = io_make_output_field("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL,
+                                 gparts, potential);
+}
+
+#endif /* SWIFT_POTENTIAL_GRAVITY_IO_H */
diff --git a/src/gravity/Potential/gravity_part.h b/src/gravity/Potential/gravity_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..252c18a4dc63c9cea4211ed8ab23eb692f064f00
--- /dev/null
+++ b/src/gravity/Potential/gravity_part.h
@@ -0,0 +1,80 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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_POTENTIAL_GRAVITY_PART_H
+#define SWIFT_POTENTIAL_GRAVITY_PART_H
+
+/* Gravity particle. */
+struct gpart {
+
+  /*! Particle ID. If negative, it is the negative offset of the #part with
+     which this gpart is linked. */
+  long long id_or_neg_offset;
+
+  /*! Particle position. */
+  double x[3];
+
+  /*! Particle velocity. */
+  float v_full[3];
+
+  /*! Particle acceleration. */
+  float a_grav[3];
+
+  /*! Particle mass. */
+  float mass;
+
+  /*! Gravitational potential */
+  float potential;
+
+  /*! Time-step length */
+  timebin_t time_bin;
+
+  /*! Type of the #gpart (DM, gas, star, ...) */
+  enum part_type type;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+  /* Numer of gparts this gpart interacted with */
+  long long num_interacted;
+
+  /* Time of the last drift */
+  integertime_t ti_drift;
+
+  /* Time of the last kick */
+  integertime_t ti_kick;
+
+#endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+  /*! Acceleration taken from the mesh only */
+  float a_grav_PM[3];
+
+  /*! Potential taken from the mesh only */
+  float potential_PM;
+
+  /* Brute-force particle acceleration. */
+  double a_grav_exact[3];
+
+  /* Brute-force particle potential. */
+  double potential_exact;
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_POTENTIAL_GRAVITY_PART_H */
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 9101fc763d69ace10c57536ebf59f655744b212b..821f044429b445c28ff8ae39b8dc65304dd2b42d 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -429,7 +429,7 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
  * @param gparts The #gpart array to write to.
  * @param gcount The number of particles to write.
  */
-__attribute__((always_inline)) INLINE void gravity_cache_write_back(
+__attribute__((always_inline)) INLINE static void gravity_cache_write_back(
     const struct gravity_cache *c, struct gpart *restrict gparts,
     const int gcount) {
 
@@ -446,7 +446,7 @@ __attribute__((always_inline)) INLINE void gravity_cache_write_back(
       gparts[i].a_grav[0] += a_x[i];
       gparts[i].a_grav[1] += a_y[i];
       gparts[i].a_grav[2] += a_z[i];
-      gparts[i].potential += pot[i];
+      gravity_add_comoving_potential(&gparts[i], pot[i]);
     }
   }
 }
diff --git a/src/gravity_iact.h b/src/gravity_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..2aaf7219d0851058cba92d4686ed989572dbcaa4
--- /dev/null
+++ b/src/gravity_iact.h
@@ -0,0 +1,39 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 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_GRAVITY_IACT_H
+#define SWIFT_GRAVITY_IACT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "inline.h"
+#include "part.h"
+
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
+#include "./gravity/Default/gravity_iact.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_iact.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
+
+#endif
diff --git a/src/gravity_io.h b/src/gravity_io.h
index 6276e50473de64abd1014bdb36a63a14e02ca8cf..752ee906081d706865e62bae7c8f505e9ca64347 100644
--- a/src/gravity_io.h
+++ b/src/gravity_io.h
@@ -19,8 +19,19 @@
 #ifndef SWIFT_GRAVITY_IO_H
 #define SWIFT_GRAVITY_IO_H
 
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
 #include "./const.h"
 
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity_io.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_io.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
 
 #endif /* SWIFT_GRAVITY_IO_H */
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 22856c25fbbaebd1dd74c5592ce8fd914e76b61c..fc1ce1d62e02c32d44667d602448fc4eb3a65344 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -112,6 +112,8 @@ void gravity_update(struct gravity_props *p, const struct cosmology *cosmo) {
 
 void gravity_props_print(const struct gravity_props *p) {
 
+  message("Self-gravity scheme: %s", GRAVITY_IMPLEMENTATION);
+
   message("Self-gravity scheme: FMM-MM with m-poles of order %d",
           SELF_GRAVITY_MULTIPOLE_ORDER);
 
@@ -167,6 +169,7 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
                        "Maximal physical softening length (Plummer equivalent)",
                        p->epsilon_max_physical);
   io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
+  io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
   io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
   io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
   io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c
index 49e66763f7f9ef47ca79d0024fa65b9ce7803fc1..2359b8a9cdf785bce719a1d0379d177d00328b9e 100644
--- a/src/mesh_gravity.c
+++ b/src/mesh_gravity.c
@@ -254,7 +254,7 @@ void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, int N, double fac,
   /* ---- */
 
   /* Store things back */
-  gp->potential += p;
+  gravity_add_comoving_potential(gp, p);
   gp->a_grav[0] += fac * a[0];
   gp->a_grav[1] += fac * a[1];
   gp->a_grav[2] += fac * a[2];
@@ -278,13 +278,14 @@ void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, int N, double fac,
  * Note that there is no multiplication by G_newton at this stage.
  *
  * @param mesh The #pm_mesh used to store the potential.
- * @param e The #engine from which to compute the forces.
+ * @param s The #space containing the particles.
+ * @param verbose Are we talkative?
  */
-void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
+void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct space* s,
+                               int verbose) {
 
 #ifdef HAVE_FFTW
 
-  const struct space* s = e->s;
   const double r_s = mesh->r_s;
   const double box_size = s->dim[0];
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
@@ -313,14 +314,24 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
   fftw_plan inverse_plan = fftw_plan_dft_c2r_3d(
       N, N, N, frho, rho, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
 
+  const ticks tic = getticks();
+
+  /* Zero everything */
   bzero(rho, N * N * N * sizeof(double));
 
   /* Do a CIC mesh assignment of the gparts */
-  for (size_t i = 0; i < e->s->nr_gparts; ++i)
-    gpart_to_mesh_CIC(&e->s->gparts[i], rho, N, cell_fac, dim);
+  for (size_t i = 0; i < s->nr_gparts; ++i)
+    gpart_to_mesh_CIC(&s->gparts[i], rho, N, cell_fac, dim);
+
+  if (verbose)
+    message("gpart assignment took %.3f %s.",
+            clocks_from_ticks(getticks() - tic), clocks_getunit());
 
+  /* message("\n\n\n DENSITY"); */
   /* print_array(rho, N); */
 
+  const ticks tic2 = getticks();
+
   /* Fourier transform to go to magic-land */
   fftw_execute(forward_plan);
 
@@ -396,6 +407,10 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
   /* Let's store it in the structure */
   mesh->potential = rho;
 
+  if (verbose)
+    message("Fourier-space PM took %.3f %s.",
+            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+
   /* message("\n\n\n POTENTIAL"); */
   /* print_array(potential, N); */
 
diff --git a/src/mesh_gravity.h b/src/mesh_gravity.h
index 8cc4cd753ccfb24f80f997dd389f2e5569c7dce1..c512a53ca349816caf4c666c6f504dd4b717bcb7 100644
--- a/src/mesh_gravity.h
+++ b/src/mesh_gravity.h
@@ -27,7 +27,7 @@
 #include "restart.h"
 
 /* Forward declarations */
-struct engine;
+struct space;
 struct gpart;
 
 /**
@@ -66,7 +66,8 @@ struct pm_mesh {
 void pm_mesh_init(struct pm_mesh *mesh, const struct gravity_props *props,
                   double dim[3]);
 void pm_mesh_init_no_mesh(struct pm_mesh *mesh, double dim[3]);
-void pm_mesh_compute_potential(struct pm_mesh *mesh, const struct engine *e);
+void pm_mesh_compute_potential(struct pm_mesh *mesh, const struct space *s,
+                               int verbose);
 void pm_mesh_interpolate_forces(const struct pm_mesh *mesh,
                                 const struct engine *e, struct gpart *gparts,
                                 int gcount);
diff --git a/src/multipole.h b/src/multipole.h
index ac746728ee39d597c92bba8c260140c10adbddd1..e0e6da32a2950d7fce164b2abc422302b7c7de5e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -31,6 +31,7 @@
 #include "align.h"
 #include "const.h"
 #include "error.h"
+#include "gravity.h"
 #include "gravity_derivatives.h"
 #include "gravity_properties.h"
 #include "gravity_softened_derivatives.h"
@@ -107,10 +108,16 @@ struct grav_tensor {
 
 struct multipole {
 
-  /* Bulk velocity */
+  /*! Bulk velocity */
   float vel[3];
 
-  /* 0th order terms */
+  /*! Maximal velocity along each axis of all #gpart */
+  float max_delta_vel[3];
+
+  /*! Minimal velocity along each axis of all #gpart */
+  float min_delta_vel[3];
+
+  /* 0th order term */
   float M_000;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -214,14 +221,15 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
 /**
  * @brief Drifts a #multipole forward in time.
  *
+ * This uses a first-order approximation in time. We only move the CoM
+ * using the bulk velocity measured at the last rebuild.
+ *
  * @param m The #multipole.
  * @param dt The drift time-step.
- * @param x_diff The maximal distance moved by any particle since the last
- * rebuild.
  */
-INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
-                                 float x_diff) {
+INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
 
+  /* Motion of the centre of mass */
   const double dx = m->m_pole.vel[0] * dt;
   const double dy = m->m_pole.vel[1] * dt;
   const double dz = m->m_pole.vel[2] * dt;
@@ -231,8 +239,36 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
   m->CoM[1] += dy;
   m->CoM[2] += dz;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (m->m_pole.vel[0] > m->m_pole.max_delta_vel[0])
+    error("Invalid maximal velocity");
+  if (m->m_pole.vel[0] < m->m_pole.min_delta_vel[0])
+    error("Invalid minimal velocity");
+  if (m->m_pole.vel[1] > m->m_pole.max_delta_vel[1])
+    error("Invalid maximal velocity");
+  if (m->m_pole.vel[1] < m->m_pole.min_delta_vel[1])
+    error("Invalid minimal velocity");
+  if (m->m_pole.vel[2] > m->m_pole.max_delta_vel[2])
+    error("Invalid maximal velocity");
+  if (m->m_pole.vel[2] < m->m_pole.min_delta_vel[2])
+    error("Invalid minimal velocity");
+#endif
+
+  /* Maximal distance covered by any particle */
+  float dv[3];
+  dv[0] = max(m->m_pole.max_delta_vel[0] - m->m_pole.vel[0],
+              m->m_pole.vel[0] - m->m_pole.min_delta_vel[0]);
+  dv[1] = max(m->m_pole.max_delta_vel[1] - m->m_pole.vel[1],
+              m->m_pole.vel[1] - m->m_pole.min_delta_vel[1]);
+  dv[2] = max(m->m_pole.max_delta_vel[2] - m->m_pole.vel[2],
+              m->m_pole.vel[2] - m->m_pole.min_delta_vel[2]);
+
+  const float max_delta_vel =
+      sqrt(dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]);
+  const float x_diff = max_delta_vel * dt;
+
   /* Conservative change in maximal radius containing all gpart */
-  m->r_max = m->r_max_rebuild + x_diff;
+  m->r_max += x_diff;
 }
 
 /**
@@ -469,16 +505,8 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
 INLINE static void gravity_multipole_add(struct multipole *ma,
                                          const struct multipole *mb) {
 
-  const float M_000 = ma->M_000 + mb->M_000;
-  const float inv_M_000 = 1.f / M_000;
-
-  /* Add the bulk velocities */
-  ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
-  ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
-  ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;
-
-  /* Add 0th order terms */
-  ma->M_000 = M_000;
+  /* Add 0th order term */
+  ma->M_000 += mb->M_000;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   /* Add 1st order terms */
@@ -554,11 +582,6 @@ INLINE static void gravity_multipole_add(struct multipole *ma,
 #error "Missing implementation for order >5"
 #endif
 
-  // MATTHIEU
-  ma->M_100 = 0.f;
-  ma->M_010 = 0.f;
-  ma->M_001 = 0.f;
-
 #ifdef SWIFT_DEBUG_CHECKS
   ma->num_gpart += mb->num_gpart;
 #endif
@@ -1025,6 +1048,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
 
   /* Prepare some local counters */
   double r_max2 = 0.;
+  float max_delta_vel[3] = {0., 0., 0.};
+  float min_delta_vel[3] = {0., 0., 0.};
+
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   double M_100 = 0., M_010 = 0., M_001 = 0.;
 #endif
@@ -1067,6 +1093,16 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
     /* Maximal distance CoM<->gpart */
     r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
 
+    /* Store the vector of the maximal vel difference */
+    max_delta_vel[0] = max(gparts[k].v_full[0], max_delta_vel[0]);
+    max_delta_vel[1] = max(gparts[k].v_full[1], max_delta_vel[1]);
+    max_delta_vel[2] = max(gparts[k].v_full[2], max_delta_vel[2]);
+
+    /* Store the vector of the minimal vel difference */
+    min_delta_vel[0] = min(gparts[k].v_full[0], min_delta_vel[0]);
+    min_delta_vel[1] = min(gparts[k].v_full[1], min_delta_vel[1]);
+    min_delta_vel[2] = min(gparts[k].v_full[2], min_delta_vel[2]);
+
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
     const double m = gparts[k].mass;
 
@@ -1149,7 +1185,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
   }
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-  M_100 = M_010 = M_001 = 0.f; /* Matthieu */
+
+  /* We know the first-order multipole (dipole) is 0. */
+  M_100 = M_010 = M_001 = 0.f;
 #endif
 
   /* Store the data on the multipole. */
@@ -1161,6 +1199,12 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
   multi->m_pole.vel[0] = vel[0];
   multi->m_pole.vel[1] = vel[1];
   multi->m_pole.vel[2] = vel[2];
+  multi->m_pole.max_delta_vel[0] = max_delta_vel[0];
+  multi->m_pole.max_delta_vel[1] = max_delta_vel[1];
+  multi->m_pole.max_delta_vel[2] = max_delta_vel[2];
+  multi->m_pole.min_delta_vel[0] = min_delta_vel[0];
+  multi->m_pole.min_delta_vel[1] = min_delta_vel[1];
+  multi->m_pole.min_delta_vel[2] = min_delta_vel[2];
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
@@ -1259,10 +1303,6 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
 INLINE static void gravity_M2M(struct multipole *m_a,
                                const struct multipole *m_b,
                                const double pos_a[3], const double pos_b[3]) {
-  /* Shift bulk velocity */
-  m_a->vel[0] = m_b->vel[0];
-  m_a->vel[1] = m_b->vel[1];
-  m_a->vel[2] = m_b->vel[2];
 
   /* Shift 0th order term */
   m_a->M_000 = m_b->M_000;
@@ -2380,7 +2420,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
   gp->a_grav[0] += a_grav[0];
   gp->a_grav[1] += a_grav[1];
   gp->a_grav[2] += a_grav[2];
-  gp->potential += pot;
+  gravity_add_comoving_potential(gp, pot);
 }
 
 /**
diff --git a/src/part.h b/src/part.h
index 03ec331cb17b95b0133be568d6e857a44d1eaf73..bca84cc0212e79e15ffbeeeb0bbcfc714d5481be 100644
--- a/src/part.h
+++ b/src/part.h
@@ -77,7 +77,13 @@
 #endif
 
 /* Import the right gravity particle definition */
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity_part.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_part.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
 
 /* Import the right star particle definition */
 #include "./stars/Default/star_part.h"
diff --git a/src/physical_constants_cgs.h b/src/physical_constants_cgs.h
index 3eeb664ab095fddf1d3451ddd3451a28fc85bd9b..40eef2c992e819e01980cbcbd7ea7f05721e93cf 100644
--- a/src/physical_constants_cgs.h
+++ b/src/physical_constants_cgs.h
@@ -93,6 +93,6 @@ const double const_solar_mass_cgs = 1.98848e33;
 const double const_earth_mass_cgs = 5.9724e27;
 
 /*! Temperature of the CMB at present day [K] */
-const double const_T_CMB_0_cgs = 2.72556;
+const double const_T_CMB_0_cgs = 2.7255;
 
 #endif /* SWIFT_PHYSICAL_CONSTANTS_CGS_H */
diff --git a/src/runner.c b/src/runner.c
index fe8e7483e24406bf23bfec6a684814c664158a69..c444fe9e1796f8a03fc44e2be155a6dcb364e63d 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1663,11 +1663,11 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
   struct gpart *restrict gparts = c->gparts;
   struct spart *restrict sparts = c->sparts;
   const int periodic = s->periodic;
-  const int with_cosmology = e->policy & engine_policy_cosmology;
-  const float Omega_m = e->cosmology->Omega_m;
-  const float H0 = e->cosmology->H0;
   const float G_newton = e->physical_constants->const_newton_G;
-  const float rho_crit0 = 3.f * H0 * H0 / (8.f * M_PI * G_newton);
+  const double r_s = e->mesh->r_s;
+  const double volume = s->dim[0] * s->dim[1] * s->dim[2];
+  const float potential_normalisation =
+      4. * M_PI * e->total_mass * r_s * r_s / volume;
 
   TIMER_TIC;
 
@@ -1702,18 +1702,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
       if (gpart_is_active(gp, e)) {
 
         /* Finish the force calculation */
-        gravity_end_force(gp, G_newton);
-
-        /* Apply periodic BC contribution to the potential */
-        if (with_cosmology && periodic) {
-          const float mass = gravity_get_mass(gp);
-          const float mass2 = mass * mass;
-
-          /* This correction term matches the one used in Gadget-2 */
-          /* The numerical constant is taken from Hernquist, Bouchet & Suto 1991
-           */
-          gp->potential -= 2.8372975f * cbrtf(mass2 * Omega_m * rho_crit0);
-        }
+        gravity_end_force(gp, G_newton, potential_normalisation, periodic);
 
 #ifdef SWIFT_NO_GRAVITY_BELOW_ID
 
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0dac3ce940d4d09c4e1f6a546a379a7b9174a1b1..142601a2181432b039fec8242471dcb4edc5f93a 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -25,6 +25,7 @@
 #include "cell.h"
 #include "gravity.h"
 #include "gravity_cache.h"
+#include "gravity_iact.h"
 #include "inline.h"
 #include "part.h"
 #include "space_getsid.h"
@@ -1010,15 +1011,6 @@ static INLINE void runner_doself_grav_pp_truncated(
       runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
                                     r_s_inv, &f_ij, &pot_ij);
 
-      /* if (e->s->parts[-gparts[pid].id_or_neg_offset].id == ICHECK) { */
-      /*   if (pjd < gcount) */
-      /*     message("Interacting with particle ID= %lld f_ij=%e", */
-      /*             e->s->parts[-gparts[pjd].id_or_neg_offset].id, f_ij); */
-      /*   // else */
-      /*   //  message("Interacting with particle ID= %lld (padded) f_ij=%e", */
-      /*   //  e->s->parts[-gparts[pjd].id_or_neg_offset].id, f_ij); */
-      /* } */
-
       /* Store it back */
       a_x += f_ij * dx;
       a_y += f_ij * dy;
diff --git a/src/space.c b/src/space.c
index 20d2159f015790542c7615302058d17d71fbf70b..6f98e788e9625c1cc872f59c58a8bf87b7b2cfa8 100644
--- a/src/space.c
+++ b/src/space.c
@@ -169,7 +169,6 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->force = NULL;
     c->grav = NULL;
     c->dx_max_part = 0.0f;
-    c->dx_max_gpart = 0.0f;
     c->dx_max_sort = 0.0f;
     c->sorted = 0;
     c->count = 0;
@@ -1732,6 +1731,7 @@ void space_split_recursive(struct space *s, struct cell *c,
   const int count = c->count;
   const int gcount = c->gcount;
   const int scount = c->scount;
+  const int with_gravity = s->gravity;
   const int depth = c->depth;
   int maxdepth = 0;
   float h_max = 0.0f;
@@ -1793,8 +1793,9 @@ void space_split_recursive(struct space *s, struct cell *c,
   }
 
   /* Split or let it be? */
-  if (count > space_splitsize || gcount > space_splitsize ||
-      scount > space_splitsize) {
+  if ((with_gravity && gcount > space_splitsize) ||
+      (!with_gravity &&
+       (count > space_splitsize || scount > space_splitsize))) {
 
     /* No longer just a leaf. */
     c->split = 1;
@@ -1823,7 +1824,6 @@ void space_split_recursive(struct space *s, struct cell *c,
       cp->split = 0;
       cp->h_max = 0.f;
       cp->dx_max_part = 0.f;
-      cp->dx_max_gpart = 0.f;
       cp->dx_max_sort = 0.f;
       cp->nodeID = c->nodeID;
       cp->parent = c;
@@ -1888,22 +1888,53 @@ void space_split_recursive(struct space *s, struct cell *c,
       /* Reset everything */
       gravity_reset(c->multipole);
 
-      /* Compute CoM of all progenies */
+      /* Compute CoM and bulk velocity from all progenies */
       double CoM[3] = {0., 0., 0.};
+      double vel[3] = {0., 0., 0.};
+      float max_delta_vel[3] = {0.f, 0.f, 0.f};
+      float min_delta_vel[3] = {0.f, 0.f, 0.f};
       double mass = 0.;
 
       for (int k = 0; k < 8; ++k) {
         if (c->progeny[k] != NULL) {
           const struct gravity_tensors *m = c->progeny[k]->multipole;
+
+          mass += m->m_pole.M_000;
+
           CoM[0] += m->CoM[0] * m->m_pole.M_000;
           CoM[1] += m->CoM[1] * m->m_pole.M_000;
           CoM[2] += m->CoM[2] * m->m_pole.M_000;
-          mass += m->m_pole.M_000;
+
+          vel[0] += m->m_pole.vel[0] * m->m_pole.M_000;
+          vel[1] += m->m_pole.vel[1] * m->m_pole.M_000;
+          vel[2] += m->m_pole.vel[2] * m->m_pole.M_000;
+
+          max_delta_vel[0] = max(m->m_pole.max_delta_vel[0], max_delta_vel[0]);
+          max_delta_vel[1] = max(m->m_pole.max_delta_vel[1], max_delta_vel[1]);
+          max_delta_vel[2] = max(m->m_pole.max_delta_vel[2], max_delta_vel[2]);
+
+          min_delta_vel[0] = min(m->m_pole.min_delta_vel[0], min_delta_vel[0]);
+          min_delta_vel[1] = min(m->m_pole.min_delta_vel[1], min_delta_vel[1]);
+          min_delta_vel[2] = min(m->m_pole.min_delta_vel[2], min_delta_vel[2]);
         }
       }
-      c->multipole->CoM[0] = CoM[0] / mass;
-      c->multipole->CoM[1] = CoM[1] / mass;
-      c->multipole->CoM[2] = CoM[2] / mass;
+
+      /* Final operation on the CoM and bulk velocity */
+      const double inv_mass = 1. / mass;
+      c->multipole->CoM[0] = CoM[0] * inv_mass;
+      c->multipole->CoM[1] = CoM[1] * inv_mass;
+      c->multipole->CoM[2] = CoM[2] * inv_mass;
+      c->multipole->m_pole.vel[0] = vel[0] * inv_mass;
+      c->multipole->m_pole.vel[1] = vel[1] * inv_mass;
+      c->multipole->m_pole.vel[2] = vel[2] * inv_mass;
+
+      /* Min max velocity along each axis */
+      c->multipole->m_pole.max_delta_vel[0] = max_delta_vel[0];
+      c->multipole->m_pole.max_delta_vel[1] = max_delta_vel[1];
+      c->multipole->m_pole.max_delta_vel[2] = max_delta_vel[2];
+      c->multipole->m_pole.min_delta_vel[0] = min_delta_vel[0];
+      c->multipole->m_pole.min_delta_vel[1] = min_delta_vel[1];
+      c->multipole->m_pole.min_delta_vel[2] = min_delta_vel[2];
 
       /* Now shift progeny multipoles and add them up */
       struct multipole temp;
@@ -1925,6 +1956,7 @@ void space_split_recursive(struct space *s, struct cell *c,
           r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
         }
       }
+
       /* Alternative upper limit of max CoM<->gpart distance */
       const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
                             ? c->multipole->CoM[0] - c->loc[0]
@@ -1945,6 +1977,11 @@ void space_split_recursive(struct space *s, struct cell *c,
       c->multipole->CoM_rebuild[1] = c->multipole->CoM[1];
       c->multipole->CoM_rebuild[2] = c->multipole->CoM[2];
 
+      /* We know the first-order multipole (dipole) is 0. */
+      c->multipole->m_pole.M_100 = 0.f;
+      c->multipole->m_pole.M_010 = 0.f;
+      c->multipole->m_pole.M_001 = 0.f;
+
     } /* Deal with gravity */
   }   /* Split or let it be? */
 
@@ -1977,7 +2014,7 @@ void space_split_recursive(struct space *s, struct cell *c,
       xparts[k].x_diff[2] = 0.f;
     }
 
-    /* gparts: Get dt_min/dt_max, reset x_diff. */
+    /* gparts: Get dt_min/dt_max. */
     for (int k = 0; k < gcount; k++) {
 #ifdef SWIFT_DEBUG_CHECKS
       if (gparts[k].time_bin == time_bin_inhibited)
@@ -1985,9 +2022,6 @@ void space_split_recursive(struct space *s, struct cell *c,
 #endif
       gravity_time_bin_min = min(gravity_time_bin_min, gparts[k].time_bin);
       gravity_time_bin_max = max(gravity_time_bin_max, gparts[k].time_bin);
-      gparts[k].x_diff[0] = 0.f;
-      gparts[k].x_diff[1] = 0.f;
-      gparts[k].x_diff[2] = 0.f;
     }
 
     /* sparts: Get dt_min/dt_max */
diff --git a/theory/Cosmology/coordinates.tex b/theory/Cosmology/coordinates.tex
index bc593606026217f345e2f77d90cee3f6632b3cef..22d0da29330f850df7fc3bdb979b49c70670d96c 100644
--- a/theory/Cosmology/coordinates.tex
+++ b/theory/Cosmology/coordinates.tex
@@ -38,11 +38,14 @@ $\Psi \equiv \frac{1}{2}a\dot{a}\mathbf{r}_i^2$ and obtain
   \Lag &= \frac{1}{2}m_ia^2 \dot{\mathbf{r}}_i^2 -
   \frac{1}{\gamma-1}m_iA_i'\left(\frac{\rho_i'}{a^3}\right)^{\gamma-1}
   -\frac{\phi'}{a},\\
-  \phi' &= a\phi + \frac{1}{2}a^2\ddot{a}\mathbf{r}_i'^2.\nonumber
+  \phi' &= a\phi + \frac{1}{2}a^2\ddot{a}\mathbf{r}_i'^2,\nonumber
 \end{align}
-Finally, we introduce the velocities
-$\mathbf{v}' \equiv a^2\dot{\mathbf{r}'}$ that are used internally by
-the code. Note that these velocities \emph{do not} have a physical
+and call $\phi'$ the peculiar potential.  Finally, we introduce the
+velocities used internally by the code:
+\begin{equation}
+  \mathbf{v}' \equiv a^2\dot{\mathbf{r}'}.
+\end{equation}
+Note that these velocities \emph{do not} have a physical
 interpretation. We caution that they are not the peculiar velocities
 ($\mathbf{v}_{\rm p} \equiv a\dot{\mathbf{r}'} =
 \frac{1}{a}\mathbf{v}'$), nor the Hubble flow
@@ -59,7 +62,9 @@ codes\footnote{One additional inconvenience of our choice of
   will hence read
   $v_{\rm sig} = a\dot{\mathbf{r}'} + c = \frac{1}{a} \left(
     |\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}.
-% This choice implies that $\dot{v}' = a \ddot{r}$.
+This choice implies that $\dot{\mathbf{v}'} = 2H\mathbf{v}' +  a^2\ddot{\mathbf{r}'}$.
+
+\subsubsection{SPH equations}
 
 Using the SPH definition of density,
 $\rho_i' = \sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
@@ -134,3 +139,4 @@ then
 \end{equation}
 indicating that the entropy evolves with the same scale-factor
 dependence as the comoving positions (eq.~\ref{eq:cosmo_eom_r}).
+
diff --git a/theory/Multipoles/exact_forces.tex b/theory/Multipoles/exact_forces.tex
index 2602c2db1c06eb496d4c450e0d33ee87d28831e8..96446080b3d5cb2d3a6c9cdb2697f50e109432e6 100644
--- a/theory/Multipoles/exact_forces.tex
+++ b/theory/Multipoles/exact_forces.tex
@@ -13,8 +13,8 @@ In the case where periodic boundary conditions are used, we apply the
 the forces of all the infinite periodic replications of the particle
 distribution. We use the approximation to the infinite series of terms
 proposed by \cite{Hernquist1991}\footnote{Note, however, that there is
-a typo in their formula for the force correction terms. The correct
-expression is given by \cite{Klessen1997} \citep[see
-also][]{Hubber2011}.}, which we tabulate in one octant using 64
-equally spaced bins along each spatial direction spanning and the
-range $[0,L]$, where $L$ is the side-length of the box.
+  a typo in their formula for the force correction terms
+  (Eq. 2.14b). The correct expression is given by \cite{Klessen1997}
+  \citep[see also][]{Hubber2011}.}, which we tabulate in one octant
+using 64 equally spaced bins along each spatial direction spanning and
+the range $[0,L]$, where $L$ is the side-length of the box.
diff --git a/theory/Multipoles/mesh_summary.tex b/theory/Multipoles/mesh_summary.tex
index a4b0b3a1bbc23fb43a978ba040d25e5f29ed384f..19524ee21b9ef85e45927182d6632dbf17ab3275 100644
--- a/theory/Multipoles/mesh_summary.tex
+++ b/theory/Multipoles/mesh_summary.tex
@@ -97,3 +97,21 @@ The truncation function in Fourier space reads
 \caption{cc}
 \label{fig:fmm:potential_long}
 \end{figure}
+
+
+\subsubsection{Normalisation of the potential}
+
+The gravitational potential computed by the combination of the mesh
+and multipole method needs to be corrected to obtain the zero point of
+the peculiar potential. Learning from the Ewald summation technique
+(see Sec.~\ref{ssec:exact_forces}), we notice that of the three terms
+entering the calculation of the potential with periodic boundary
+conditions (see Eq. 2.11 of \cite{Hernquist1991}), the last two are
+computed by the mesh and tree respectively. The only reamining term is
+the constant
+\begin{equation}
+  \frac{4\pi r_s^2}{L^3}\sum_a m_a,
+\end{equation}
+where the sum runs over all particles and $L$ is the comoving
+side-length of the simulation volume. This constant is added to the
+potential of all the particles in the simulation.