diff --git a/src/multipole.c b/src/multipole.c
new file mode 100644
index 0000000000000000000000000000000000000000..6a3a52c3e74ac74e7ffa3cfd487d82b454cc9472
--- /dev/null
+++ b/src/multipole.c
@@ -0,0 +1,190 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2013 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <float.h>
+#include <limits.h>
+
+/* MPI headers. */
+#ifdef WITH_MPI
+ #include <mpi.h>
+#endif
+
+/* Local headers. */
+#include "error.h"
+#include "const.h"
+#include "cycle.h"
+#include "atomic.h"
+#include "lock.h"
+#include "space.h"
+#include "part.h"
+#include "multipole.h"
+#include "cell.h"
+
+
+/**
+ * @brief Merge two multipoles.
+ *
+ * @param ma The #multipole which will contain the merged result.
+ * @param mb The other #multipole.
+ */
+
+void multipole_merge ( struct multipole *ma , struct multipole *mb ) {
+
+ #if multipole_order == 1
+
+ /* Correct the position. */
+ float mma = ma->coeffs[0], mmb = mb->coeffs[0];
+ float w = 1.0f / ( mma + mmb );
+ for ( int k = 0 ; k < 3 ; k++ )
+ ma->x[k] = ( ma->x[k]*mma + mb->x[k]*mmb ) * w;
+
+ /* Add the particle to the moments. */
+ ma->coeffs[0] = mma + mmb;
+
+ #else
+ #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ #endif
+
+ }
+
+
+/**
+ * @brief Add a particle to the given multipole.
+ *
+ * @param m The #multipole.
+ * @param p The #gpart.
+ */
+
+void multipole_addpart ( struct multipole *m , struct gpart *p ) {
+
+ #if multipole_order == 1
+
+ /* Correct the position. */
+ float mm = m->coeffs[0], mp = p->mass;
+ float w = 1.0f / ( mm + mp );
+ for ( int k = 0 ; k < 3 ; k++ )
+ m->x[k] = ( m->x[k]*mm + p->x[k]*mp ) * w;
+
+ /* Add the particle to the moments. */
+ m->coeffs[0] = mm + mp;
+
+ #else
+ #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ #endif
+
+ }
+
+
+/**
+ * @brief Add a group of particles to the given multipole.
+ *
+ * @param m The #multipole.
+ * @param p The #gpart array.
+ * @param N Number of parts to add.
+ */
+
+void multipole_addparts ( struct multipole *m , struct gpart *p , int N ) {
+
+ #if multipole_order == 1
+
+ /* Get the combined mass and positions. */
+ double xp[3] = { 0.0 , 0.0 , 0.0 };
+ float mp = 0.0f, w;
+ for ( int k = 0 ; k < N ; k++ ) {
+ w = p[k].mass;
+ mp += w;
+ xp[0] += p[k].x[0] * w;
+ xp[1] += p[k].x[1] * w;
+ xp[2] += p[k].x[2] * w;
+ }
+
+ /* Correct the position. */
+ float mm = m->coeffs[0];
+ w = 1.0f / ( mm + mp );
+ for ( int k = 0 ; k < 3 ; k++ )
+ m->x[k] = ( m->x[k]*mm + xp[k] ) * w;
+
+ /* Add the particle to the moments. */
+ m->coeffs[0] = mm + mp;
+
+ #else
+ #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ #endif
+
+ }
+
+
+/**
+ * @brief Init a multipole from a set of particles.
+ *
+ * @param m The #multipole.
+ * @param parts The #gparts.
+ * @param N The number of particles.
+ */
+
+void multipole_init ( struct multipole *m , struct gpart *parts , int N ) {
+
+ #if multipole_order == 1
+
+ float mass = 0.0f, w;
+ double x[3] = { 0.0 , 0.0 , 0.0 };
+ int k;
+
+ /* Collect the particle data. */
+ for ( k = 0 ; k < N ; k++ ) {
+ w = parts[k].mass;
+ mass += w;
+ x[0] += parts[k].x[0] * w;
+ x[1] += parts[k].x[1] * w;
+ x[2] += parts[k].x[2] * w;
+ }
+
+ /* Store the data on the multipole. */
+ m->coeffs[0] = mass;
+ m->x[0] = x[0] / mass;
+ m->x[1] = x[1] / mass;
+ m->x[2] = x[2] / mass;
+
+ #else
+ #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ #endif
+
+ }
+
+
+/**
+ * @brief Reset the data of a #multipole.
+ *
+ * @param m The #multipole.
+ */
+
+void multipole_reset ( struct multipole *m ) {
+
+ /* Just bzero the struct. */
+ bzero( m , sizeof(struct multipole) );
+
+ }
diff --git a/src/multipole.h b/src/multipole.h
new file mode 100644
index 0000000000000000000000000000000000000000..e3372e245637dbaec6e078ac7294902f6a22a9fe
--- /dev/null
+++ b/src/multipole.h
@@ -0,0 +1,133 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2013 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/>.
+ *
+ ******************************************************************************/
+
+/* Some constants. */
+#define multipole_order 1
+
+
+/* Multipole struct. */
+struct multipole {
+
+ /* Multipole location. */
+ double x[3];
+
+ /* Acceleration on this multipole. */
+ float a[3];
+
+ /* Multipole coefficients. */
+ float coeffs[ multipole_order*multipole_order ];
+
+ };
+
+
+/* Multipole function prototypes. */
+static void multipole_iact_mm ( struct multipole *ma , struct multipole *mb , double *shift );
+void multipole_merge ( struct multipole *ma , struct multipole *mb );
+void multipole_addpart ( struct multipole *m , struct gpart *p );
+void multipole_addparts ( struct multipole *m , struct gpart *p , int N );
+void multipole_init ( struct multipole *m , struct gpart *parts , int N );
+void multipole_reset ( struct multipole *m );
+
+
+#include <math.h>
+#include "kernel.h"
+
+/**
+ * @brief Compute the pairwise interaction between two multipoles.
+ *
+ * @param ma The first #multipole.
+ * @param mb The second #multipole.
+ * @param shift The periodicity correction.
+ */
+
+__attribute__ ((always_inline)) INLINE static void multipole_iact_mm ( struct multipole *ma , struct multipole *mb , double *shift ) {
+
+ float dx[3], ir, r, r2 = 0.0f, acc;
+ int k;
+
+ /* Compute the multipole distance. */
+ for ( k = 0 ; k < 3 ; k++ ) {
+ dx[k] = ma->x[k] - mb->x[k] - shift[k];
+ r2 += dx[k]*dx[k];
+ }
+
+ /* Compute the normalized distance vector. */
+ ir = 1.0f / sqrtf( r2 );
+ r = r2 * ir;
+
+ /* Evaluate the gravity kernel. */
+ kernel_grav_eval( r , &acc );
+
+ /* Scale the acceleration. */
+ acc *= const_G * ir * ir * ir;
+
+ /* Compute the forces on both multipoles. */
+ #if multipole_order == 1
+ float mma = ma->coeffs[0], mmb = mb->coeffs[0];
+ for ( k = 0 ; k < 3 ; k++ ) {
+ ma->a[k] -= dx[k] * acc * mmb;
+ mb->a[k] += dx[k] * acc * mma;
+ }
+ #else
+ #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ #endif
+
+ }
+
+
+/**
+ * @brief Compute the interaction of a multipole on a particle.
+ *
+ * @param m The #multipole.
+ * @param p The #gpart.
+ * @param shift The periodicity correction.
+ */
+
+__attribute__ ((always_inline)) INLINE static void multipole_iact_mp ( struct multipole *m , struct gpart *p , double *shift ) {
+
+ float dx[3], ir, r, r2 = 0.0f, acc;
+ int k;
+
+ /* Compute the multipole distance. */
+ for ( k = 0 ; k < 3 ; k++ ) {
+ dx[k] = m->x[k] - p->x[k] - shift[k];
+ r2 += dx[k]*dx[k];
+ }
+
+ /* Compute the normalized distance vector. */
+ ir = 1.0f / sqrtf( r2 );
+ r = r2 * ir;
+
+ /* Evaluate the gravity kernel. */
+ kernel_grav_eval( r , &acc );
+
+ /* Scale the acceleration. */
+ acc *= const_G * ir * ir * ir * m->coeffs[0];
+
+ /* Compute the forces on both multipoles. */
+ #if multipole_order == 1
+ for ( k = 0 ; k < 3 ; k++ )
+ p->a[k] += dx[k] * acc;
+ #else
+ #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ #endif
+
+ }
+
+