diff --git a/src/cell.h b/src/cell.h
index 0f8af822bb4050a21508a6b244ca4431d4a793b3..f67b8b7455bc2d4e22b00098fd9fb6a5bcebda8a 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -64,6 +64,9 @@ struct cell {
     /* The tasks computing this cell's sorts. */
     struct task *sorts[14];
     
+    /* The ghost task to link density to interactions. */
+    struct task *ghost;
+    
     /* Number of tasks this cell is waiting for and whether it is in use. */
     int wait;
     
diff --git a/src/gadgetsmp.h b/src/gadgetsmp.h
index b10ba58d4b6473da978919ed2685c92b432cca53..a14888a70ea2cd972b0e152c70ff7f3781fb77d3 100644
--- a/src/gadgetsmp.h
+++ b/src/gadgetsmp.h
@@ -28,5 +28,6 @@
 #include "cell.h"
 #include "space.h"
 #include "queue.h"
+#include "runner_iact.h"
 #include "runner.h"
 
diff --git a/src/part.h b/src/part.h
index 8a95bf05992b8837dda93da5910e568360bc4499..07e639a1aa59166f7ddf6fd25e84738d38ada638 100644
--- a/src/part.h
+++ b/src/part.h
@@ -26,11 +26,8 @@
 /* Data of a single particle. */
 struct part {
 
-    /* Particle position. */
-    double x[3];
-    
     /* Particle cutoff radius. */
-    float r;
+    float h;
     
     /* Particle time-step. */
     float dt;
@@ -38,8 +35,34 @@ struct part {
     /* Particle ID. */
     int id;
     
-    /* Number of pairwise interactions. */
-    double count, count_dh;
+    /* Particle density. */
+    float rho;
+    
+    /* Particle position. */
+    double x[3];
+    
+    /* Particle velocity. */
+    float v[3];
+    
+    /* Particle acceleration. */
+    float a[3];
+    
+    /* Particle pressure. */
+    float P;
+    
+    /* Particle mass. */
+    float m;
+    
+    /* Particle internal energy. */
+    float u;
+    
+    /* Change in particle energy over time. */
+    float u_dt;
+    
+    /* Derivative of the density with respect to this particle's smoothing length. */
+    float rho_dh;
+    
+    /* Particle number density. */
     int icount;
     
     } __attribute__((aligned (32)));
diff --git a/src/queue.c b/src/queue.c
index f14db6a569132a6cc50134a05bb29a1f04bdff05..fca97701093dc04a7f99f1e1b5e41e044dc76566 100644
--- a/src/queue.c
+++ b/src/queue.c
@@ -185,11 +185,11 @@ struct task *queue_gettask ( struct queue *q , int blocking , int keep ) {
                 continue;
             
             /* Different criteria for different types. */
-            if ( res->type == tid_self || (res->type == tid_sub && res->cj == NULL) ) {
+            if ( res->type == task_type_self || (res->type == task_type_sub && res->cj == NULL) ) {
                 if ( res->ci->hold || cell_locktree( res->ci ) != 0 )
                     continue;
                 }
-            else if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) ) {
+            else if ( res->type == task_type_pair || (res->type == task_type_sub && res->cj != NULL) ) {
                 if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait )
                     continue;
                 if ( cell_locktree( res->ci ) != 0 )
@@ -306,7 +306,7 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int
                 continue;
                 
             /* Get the score for this task. */
-            if ( res->type == tid_self || res->type == tid_sort || ( res->type == tid_sub && res->cj == NULL ) )
+            if ( res->type == task_type_self || res->type == task_type_sort || ( res->type == task_type_sub && res->cj == NULL ) )
                 score = ( res->ci->owner == rid );
             else
                 score = ( res->ci->owner == rid ) + ( res->cj->owner == rid );
@@ -314,11 +314,11 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int
                 continue;
             
             /* Different criteria for different types. */
-            if ( res->type == tid_self || (res->type == tid_sub && res->cj == NULL) ) {
+            if ( res->type == task_type_self || (res->type == task_type_sub && res->cj == NULL) ) {
                 if ( res->ci->hold || cell_locktree( res->ci ) != 0 )
                     continue;
                 }
-            else if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) ) {
+            else if ( res->type == task_type_pair || (res->type == task_type_sub && res->cj != NULL) ) {
                 if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait )
                     continue;
                 if ( cell_locktree( res->ci ) != 0 )
@@ -332,9 +332,9 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int
             /* If we owned a previous task, unlock it. */
             if ( ind_best >= 0 ) {
                 res = &qtasks[ qtid[ ind_best ] ];
-                if ( res->type == tid_self || res->type == tid_pair || res->type == tid_sub )
+                if ( res->type == task_type_self || res->type == task_type_pair || res->type == task_type_sub )
                     cell_unlocktree( res->ci );
-                if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) )
+                if ( res->type == task_type_pair || (res->type == task_type_sub && res->cj != NULL) )
                     cell_unlocktree( res->cj );
                 }
             
@@ -438,19 +438,19 @@ void queue_sort ( struct queue *q ) {
     for ( k = 0 ; k < q->count ; k++ ) {
         t = &q->tasks[ q->tid[k] ];
         switch ( t->type ) {
-            case tid_self:
+            case task_type_self:
                 wait[k] = t->rank;
                 weight[k] = 0; // t->ci->count * t->ci->count;
                 break;
-            case tid_pair:
+            case task_type_pair:
                 wait[k] = t->rank;
                 weight[k] = 0; // t->ci->count * t->cj->count;
                 break;
-            case tid_sub:
+            case task_type_sub:
                 wait[k] = t->rank;
                 weight[k] = 0; // (t->cj == NULL) ? t->ci->count * t->ci->count : t->ci->count * t->cj->count;
                 break;
-            case tid_sort:
+            case task_type_sort:
                 wait[k] = t->rank;
                 weight[k] = 0; // t->ci->count;
                 break;
diff --git a/src/runner.c b/src/runner.c
index f224eef17b8dc452020645c8df0fb41d06b7c393..b7f10d76e73f853a7e027e4f5cee09b56b033958 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -40,6 +40,7 @@
 #include "space.h"
 #include "queue.h"
 #include "runner.h"
+#include "runner_iact.h"
 
 /* Error macro. */
 #define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
@@ -74,6 +75,15 @@ const char runner_flip[27] = { 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1
                                0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 }; 
 
 
+/* Import the density functions. */
+#define FUNCTION density
+#include "runner_doiact.h"
+
+#undef FUNCTION
+#define FUNCTION force
+#include "runner_doiact.h"
+
+
 /** 
  * @brief Sort the tasks in topological order over all queues.
  *
@@ -115,7 +125,7 @@ void runner_ranktasks ( struct runner *r ) {
             t->rank = rank;
             s->tasks_ind[i] = t - s->tasks;
             /* printf( "runner_ranktasks: task %i of type %s has rank %i.\n" , i , 
-                (t->type == tid_self) ? "self" : (t->type == tid_pair) ? "pair" : "sort" , rank ); */
+                (t->type == task_type_self) ? "self" : (t->type == task_type_pair) ? "pair" : "sort" , rank ); */
             for ( k = 0 ; k < t->nr_unlock_tasks ; k++ )
                 t->unlock_tasks[k]->wait -= 1;
             }
@@ -129,505 +139,6 @@ void runner_ranktasks ( struct runner *r ) {
     free(tid);
     
     }
-    
-
-/**
- * @brief Compute the interactions between a cell pair.
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- */
- 
-void runner_dopair_naive ( struct runner_thread *rt , struct cell *ci , struct cell *cj ) {
-
-    struct runner *r = rt->r;
-    int pid, pjd, k, count_i = ci->count, count_j = cj->count;
-    double shift[3] = { 0.0 , 0.0 , 0.0 };
-    struct part *pi, *pj, *parts_i = ci->parts, *parts_j = cj->parts;
-    double dx[3], pix[3], ri, ri2, r2;
-    TIMER_TIC
-    
-    /* Get the relative distance between the pairs, wrapping. */
-    for ( k = 0 ; k < 3 ; k++ ) {
-        if ( cj->loc[k] - ci->loc[k] < -r->s->dim[k]/2 )
-            shift[k] = r->s->dim[k];
-        else if ( cj->loc[k] - ci->loc[k] > r->s->dim[k]/2 )
-            shift[k] = -r->s->dim[k];
-        }
-        
-    /* printf( "runner_dopair_naive: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts and shift = [ %g %g %g ].\n" ,
-        ci->loc[0] , ci->loc[1] , ci->loc[2] , cj->loc[0] , cj->loc[1] , cj->loc[2] ,
-        ci->count , cj->count , shift[0] , shift[1] , shift[2] ); fflush(stdout);
-    tic = getticks(); */
-    
-    /* Loop over the parts in ci. */
-    for ( pid = 0 ; pid < count_i ; pid++ ) {
-    
-        /* Get a hold of the ith part in ci. */
-        pi = &parts_i[ pid ];
-        for ( k = 0 ; k < 3 ; k++ )
-            pix[k] = pi->x[k] - shift[k];
-        ri = pi->r;
-        ri2 = ri * ri;
-        
-        /* Loop over the parts in cj. */
-        for ( pjd = 0 ; pjd < count_j ; pjd++ ) {
-        
-            /* Get a pointer to the jth particle. */
-            pj = &parts_j[ pjd ];
-        
-            /* Compute the pairwise distance. */
-            r2 = 0.0;
-            for ( k = 0 ; k < 3 ; k++ ) {
-                dx[k] = pix[k] - pj->x[k];
-                r2 += dx[k]*dx[k];
-                }
-                
-            /* Hit or miss? */
-            if ( r2 < ri2 || r2 < pj->r*pj->r ) {
-            
-                iact( r2 , ri , pj->r , pi , pj );
-            
-                }
-        
-            } /* loop over the parts in cj. */
-    
-        } /* loop over the parts in ci. */
-        
-    #ifdef TIMER_VERBOSE
-        printf( "runner_dopair_naive[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f) took %.3f ms.\n" , rt->id , count_i , count_j , ci->depth , ci->r_max , cj->r_max , ((double)TIMER_TOC(runner_timer_dopair)) / CPU_TPS * 1000 );
-    #else
-        TIMER_TOC(runner_timer_dopair);
-    #endif
-
-
-    }
-
-
-/**
- * @brief Compute the interactions between a cell pair.
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- */
- 
-void runner_dopair ( struct runner_thread *rt , struct cell *ci , struct cell *cj ) {
-
-    struct runner *r = rt->r;
-    int pid, pjd, k, sid;
-    double rshift, shift[3] = { 0.0 , 0.0 , 0.0 };
-    struct cell *temp;
-    struct entry *sort_i, *sort_j;
-    struct part *pi, *pj, *parts_i, *parts_j;
-    double dx[3], pix[3], pjx[3], ri, ri2, rj, rj2, r2, di, dj;
-    double ri_max, rj_max, di_max, dj_min;
-    int count_i, count_j;
-    TIMER_TIC
-    
-    /* Get the relative distance between the pairs, wrapping. */
-    for ( k = 0 ; k < 3 ; k++ ) {
-        if ( cj->loc[k] - ci->loc[k] < -r->s->dim[k]/2 )
-            shift[k] = r->s->dim[k];
-        else if ( cj->loc[k] - ci->loc[k] > r->s->dim[k]/2 )
-            shift[k] = -r->s->dim[k];
-        }
-        
-    /* Get the sorting index. */
-    for ( sid = 0 , k = 0 ; k < 3 ; k++ )
-        sid = 3*sid + ( (cj->loc[k] - ci->loc[k] + shift[k] < 0) ? 0 : (cj->loc[k] - ci->loc[k] + shift[k] > 0) ? 2 : 1 );
-
-    /* Switch the cells around? */
-    if ( runner_flip[sid] ) {
-        temp = ci; ci = cj; cj = temp;
-        for ( k = 0 ; k < 3 ; k++ )
-            shift[k] = -shift[k];
-        }
-    sid = sortlistID[sid];
-    
-    /* Get the cutoff shift. */
-    for ( rshift = 0.0 , k = 0 ; k < 3 ; k++ )
-        rshift += shift[k]*runner_shift[ 3*sid + k ];
-        
-    /* printf( "runner_dopair: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts and shift = [ %g %g %g ].\n" ,
-        ci->loc[0] , ci->loc[1] , ci->loc[2] , cj->loc[0] , cj->loc[1] , cj->loc[2] ,
-        ci->count , cj->count , shift[0] , shift[1] , shift[2] ); fflush(stdout); */
-    /* for ( ri = 0 , k = 0 ; k < ci->count ; k++ )
-        ri += ci->parts[k].r;
-    for ( rj = 0 , k = 0 ; k < cj->count ; k++ )
-        rj += cj->parts[k].r;
-    printf( "runner_dopair: avg. radii %g/%g for h=%g at depth=%i.\n" , ri/ci->count , rj/cj->count , ci->h[0] , ci->depth ); fflush(stdout); */
-    
-    /* Pick-out the sorted lists. */
-    sort_i = &ci->sort[ sid*(ci->count + 1) ];
-    sort_j = &cj->sort[ sid*(cj->count + 1) ];
-    
-    /* Get some other useful values. */
-    ri_max = ci->r_max - rshift; rj_max = cj->r_max - rshift;
-    count_i = ci->count; count_j = cj->count;
-    parts_i = ci->parts; parts_j = cj->parts;
-    di_max = sort_i[count_i-1].d - rshift;
-    dj_min = sort_j[0].d;
-    
-    /* if ( ci->split && cj->split && sid == 4 )
-        printf( "boing!\n" ); */
-    
-    /* Loop over the parts in ci. */
-    for ( pid = count_i-1 ; pid >= 0 && sort_i[pid].d + ri_max > dj_min ; pid-- ) {
-    
-        /* Get a hold of the ith part in ci. */
-        pi = &parts_i[ sort_i[ pid ].i ];
-        ri = pi->r;
-        di = sort_i[pid].d + ri - rshift;
-        if ( di < dj_min )
-            continue;
-            
-        ri2 = pi->r * pi->r;
-        for ( k = 0 ; k < 3 ; k++ )
-            pix[k] = pi->x[k] - shift[k];
-        
-        /* Loop over the parts in cj. */
-        for ( pjd = 0 ; pjd < count_j && sort_j[pjd].d < di ; pjd++ ) {
-        
-            /* Get a pointer to the jth particle. */
-            pj = &parts_j[ sort_j[pjd].i ];
-        
-            /* Compute the pairwise distance. */
-            r2 = 0.0;
-            for ( k = 0 ; k < 3 ; k++ ) {
-                dx[k] = pix[k] - pj->x[k];
-                r2 += dx[k]*dx[k];
-                }
-                
-            /* Hit or miss? */
-            if ( r2 < ri2 ) {
-            
-                iact( r2 , ri , pj->r , pi , pj );
-            
-                }
-        
-            } /* loop over the parts in cj. */
-    
-        } /* loop over the parts in ci. */
-        
-    /* printf( "runner_dopair: first half took %.3f ms...\n" , ((double)(getticks() - tic)) / CPU_TPS * 1000 );
-    tic = getticks(); */
-
-    /* Loop over the parts in cj. */
-    for ( pjd = 0 ; pjd < count_j && sort_j[pjd].d - rj_max < di_max ; pjd++ ) {
-    
-        /* Get a hold of the jth part in cj. */
-        pj = &parts_j[ sort_j[ pjd ].i ];
-        rj = pj->r;
-        dj = sort_j[pjd].d - rj - rshift;
-        if ( dj > di_max )
-            continue;
-            
-        for ( k = 0 ; k < 3 ; k++ )
-            pjx[k] = pj->x[k] + shift[k];
-        rj2 = pj->r * pj->r;
-        
-        /* Loop over the parts in ci. */
-        for ( pid = count_i-1 ; pid >= 0 && sort_i[pid].d > dj ; pid-- ) {
-        
-            /* Get a pointer to the jth particle. */
-            pi = &parts_i[ sort_i[pid].i ];
-            
-            /* Compute the pairwise distance. */
-            r2 = 0.0;
-            for ( k = 0 ; k < 3 ; k++ ) {
-                dx[k] = pi->x[k] - pjx[k];
-                r2 += dx[k]*dx[k];
-                }
-                
-            /* Hit or miss? */
-            if ( r2 < rj2 && r2 > pi->r*pi->r ) {
-            
-                iact( r2 , pi->r , rj , pi , pj );
-            
-                }
-        
-            } /* loop over the parts in cj. */
-    
-        } /* loop over the parts in ci. */
-
-    #ifdef TIMER_VERBOSE
-        printf( "runner_dopair[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f, h=%.3f) took %.3f ms.\n" , rt->id , count_i , count_j , ci->depth , ci->r_max , cj->r_max , fmax(ci->h[0],fmax(ci->h[1],ci->h[2])) , ((double)(TIMER_TOC(runner_timer_dopair))) / CPU_TPS * 1000 );
-    #else
-        TIMER_TOC(runner_timer_dopair);
-    #endif
-
-    }
-
-
-/**
- * @brief Compute the cell self-interaction.
- *
- * @param r The #runner.
- * @param c The #cell.
- */
-
-void runner_doself ( struct runner_thread *rt , struct cell *c ) {
-
-    int k, pid, pjd, count = c->count;
-    double pix[3], dx[3], ri, ri2, r2;
-    struct part *pi, *pj, *parts = c->parts;
-    TIMER_TIC
-    
-    if ( c->split )
-        error( "Split cell should not have self-interactions." );
-    
-    /* Loop over the particles in the cell. */
-    for ( pid = 0 ; pid < count ; pid++ ) {
-    
-        /* Get a pointer to the ith particle. */
-        pi = &parts[pid];
-    
-        /* Get the particle position and radius. */
-        for ( k = 0 ; k < 3 ; k++ )
-            pix[k] = pi->x[k];
-        ri = pi->r;
-        ri2 = ri * ri;
-            
-        /* Loop over the other particles .*/
-        for ( pjd = pid+1 ; pjd < count ; pjd++ ) {
-        
-            /* Get a pointer to the jth particle. */
-            pj = &parts[pjd];
-        
-            /* Compute the pairwise distance. */
-            r2 = 0.0;
-            for ( k = 0 ; k < 3 ; k++ ) {
-                dx[k] = pix[k] - pj->x[k];
-                r2 += dx[k]*dx[k];
-                }
-                
-            /* Hit or miss? */
-            if ( r2 < ri2 || r2 < pj->r*pj->r ) {
-            
-                iact( r2 , ri , pj->r , pi , pj );
-            
-                }
-        
-            } /* loop over all other particles. */
-    
-        } /* loop over all particles. */
-
-    #ifdef TIMER_VERBOSE
-        printf( "runner_doself[%02i]: %i parts at depth %i took %.3f ms.\n" , rt->id , count , c->depth , ((double)TIMER_TOC(runner_timer_doself)) / CPU_TPS * 1000 );
-    #else
-        TIMER_TOC(runner_timer_doself);
-    #endif
-
-    }
-
-
-/**
- * @brief Compute grouped sub-cell interactions
- *
- * @param r The #runner.
- * @param c The #cell.
- */
-
-void runner_dosub ( struct runner_thread *rt , struct cell *ci , struct cell *cj , int flags ) {
-
-    int j, k;
-
-    TIMER_TIC
-    
-    /* Different types of flags. */
-    switch ( flags ) {
-    
-        /* Regular sub-cell interactions of a single cell. */
-        case 0:
-            for ( j = 0 ; j < 7 ; j++ )
-                for ( k = j + 1 ; k < 8 ; k++ )
-                    if ( ci->progeny[j] != NULL && ci->progeny[k] != NULL )
-                        runner_dopair( rt , ci->progeny[j] , ci->progeny[k] );
-            break;
-            
-        case 1: /* (  1 ,  1 ,  0 ) */
-            if ( ci->progeny[6] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[0] );
-            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[1] );
-            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[0] );
-            if ( ci->progeny[7] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[1] );
-            break;
-    
-        case 3: /* (  1 ,  0 ,  1 ) */
-            if ( ci->progeny[5] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[0] );
-            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[2] );
-            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[0] );
-            if ( ci->progeny[7] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[2] );
-            break;
-                    
-        case 4: /* (  1 ,  0 ,  0 ) */
-            if ( ci->progeny[4] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[0] );
-            if ( ci->progeny[4] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[1] );
-            if ( ci->progeny[4] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[2] );
-            if ( ci->progeny[4] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[3] );
-            if ( ci->progeny[5] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[0] );
-            if ( ci->progeny[5] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[1] );
-            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[2] );
-            if ( ci->progeny[5] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[3] );
-            if ( ci->progeny[6] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[0] );
-            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[1] );
-            if ( ci->progeny[6] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[2] );
-            if ( ci->progeny[6] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[3] );
-            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[0] );
-            if ( ci->progeny[7] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[1] );
-            if ( ci->progeny[7] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[2] );
-            if ( ci->progeny[7] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[3] );
-            break;
-            
-        case 5: /* (  1 ,  0 , -1 ) */
-            if ( ci->progeny[4] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[1] );
-            if ( ci->progeny[4] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[3] );
-            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[1] );
-            if ( ci->progeny[6] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[3] );
-            break;
-                    
-        case 7: /* (  1 , -1 ,  0 ) */
-            if ( ci->progeny[4] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[2] );
-            if ( ci->progeny[4] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[4] , cj->progeny[3] );
-            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[2] );
-            if ( ci->progeny[5] != NULL && cj->progeny[3] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[3] );
-            break;
-                    
-        case 9: /* (  0 ,  1 ,  1 ) */
-            if ( ci->progeny[3] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[0] );
-            if ( ci->progeny[3] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[4] );
-            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[0] );
-            if ( ci->progeny[7] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[4] );
-            break;
-                    
-        case 10: /* (  0 ,  1 ,  0 ) */
-            if ( ci->progeny[2] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[2] , cj->progeny[0] );
-            if ( ci->progeny[2] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[2] , cj->progeny[1] );
-            if ( ci->progeny[2] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[2] , cj->progeny[4] );
-            if ( ci->progeny[2] != NULL && cj->progeny[5] != NULL )
-                runner_dopair( rt , ci->progeny[2] , cj->progeny[5] );
-            if ( ci->progeny[3] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[0] );
-            if ( ci->progeny[3] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[1] );
-            if ( ci->progeny[3] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[4] );
-            if ( ci->progeny[3] != NULL && cj->progeny[5] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[5] );
-            if ( ci->progeny[6] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[0] );
-            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[1] );
-            if ( ci->progeny[6] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[4] );
-            if ( ci->progeny[6] != NULL && cj->progeny[5] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[5] );
-            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[0] );
-            if ( ci->progeny[7] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[1] );
-            if ( ci->progeny[7] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[4] );
-            if ( ci->progeny[7] != NULL && cj->progeny[5] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[5] );
-            break;
-                    
-        case 11: /* (  0 ,  1 , -1 ) */
-            if ( ci->progeny[2] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[2] , cj->progeny[1] );
-            if ( ci->progeny[2] != NULL && cj->progeny[5] != NULL )
-                runner_dopair( rt , ci->progeny[2] , cj->progeny[5] );
-            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[1] );
-            if ( ci->progeny[6] != NULL && cj->progeny[5] != NULL )
-                runner_dopair( rt , ci->progeny[6] , cj->progeny[5] );
-            break;
-                    
-        case 12: /* (  0 ,  0 ,  1 ) */
-            if ( ci->progeny[1] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[1] , cj->progeny[0] );
-            if ( ci->progeny[1] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[1] , cj->progeny[2] );
-            if ( ci->progeny[1] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[1] , cj->progeny[4] );
-            if ( ci->progeny[1] != NULL && cj->progeny[6] != NULL )
-                runner_dopair( rt , ci->progeny[1] , cj->progeny[6] );
-            if ( ci->progeny[3] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[0] );
-            if ( ci->progeny[3] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[2] );
-            if ( ci->progeny[3] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[4] );
-            if ( ci->progeny[3] != NULL && cj->progeny[6] != NULL )
-                runner_dopair( rt , ci->progeny[3] , cj->progeny[6] );
-            if ( ci->progeny[5] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[0] );
-            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[2] );
-            if ( ci->progeny[5] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[4] );
-            if ( ci->progeny[5] != NULL && cj->progeny[6] != NULL )
-                runner_dopair( rt , ci->progeny[5] , cj->progeny[6] );
-            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[0] );
-            if ( ci->progeny[7] != NULL && cj->progeny[2] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[2] );
-            if ( ci->progeny[7] != NULL && cj->progeny[4] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[4] );
-            if ( ci->progeny[7] != NULL && cj->progeny[6] != NULL )
-                runner_dopair( rt , ci->progeny[7] , cj->progeny[6] );
-            break;
-                
-        }
-    
-
-    #ifdef TIMER_VERBOSE
-        printf( "runner_dosub[%02i]: flags=%i at depth %i took %.3f ms.\n" , rt->id , flags , ci->depth , ((double)TIMER_TOC(runner_timer_dosub)) / CPU_TPS * 1000 );
-    #else
-        TIMER_TOC(runner_timer_dosub);
-    #endif
-
-    }
 
 
 /**
@@ -1132,24 +643,41 @@ void *runner_main ( void *data ) {
             
             /* Different types of tasks... */
             switch ( t->type ) {
-                case tid_self:
-                    runner_doself( rt , ci );
+                case task_type_self:
+                    if ( t->subtype == task_subtype_density )
+                        runner_doself_density( rt , ci );
+                    else if ( t->subtype == task_subtype_force )
+                        runner_doself_force( rt , ci );
+                    else
+                        error( "Unknown task subtype." );
                     cell_unlocktree( ci );
                     break;
-                case tid_pair:
-                    runner_dopair( rt , ci , cj );
+                case task_type_pair:
+                    if ( t->subtype == task_subtype_density )
+                        runner_dopair_density( rt , ci , cj );
+                    else if ( t->subtype == task_subtype_force )
+                        runner_dopair_force( rt , ci , cj );
+                    else
+                        error( "Unknown task subtype." );
                     cell_unlocktree( ci );
                     cell_unlocktree( cj );
                     break;
-                case tid_sort:
+                case task_type_sort:
                     runner_dosort( rt , ci , t->flags );
                     break;
-                case tid_sub:
-                    runner_dosub( rt , ci , cj , t->flags );
+                case task_type_sub:
+                    if ( t->subtype == task_subtype_density )
+                        runner_dosub_density( rt , ci , cj , t->flags );
+                    else if ( t->subtype == task_subtype_force )
+                        runner_dosub_force( rt , ci , cj , t->flags );
+                    else
+                        error( "Unknown task subtype." );
                     cell_unlocktree( ci );
                     if ( cj != NULL )
                         cell_unlocktree( cj );
                     break;
+                case task_type_ghost:
+                    break;
                 default:
                     error( "Unknown task type." );
                 }
@@ -1283,7 +811,7 @@ void runner_init ( struct runner *r , struct space *s , int nr_threads , int nr_
         
     /* Fill the queues (round-robin). */
     for ( k = 0 ; k < s->nr_tasks ; k++ ) {
-        if ( s->tasks[ s->tasks_ind[k] ].type == tid_none)
+        if ( s->tasks[ s->tasks_ind[k] ].type == task_type_none )
             continue;
         // qid = 0;
         // qid = k % nrq;
diff --git a/src/runner.h b/src/runner.h
index c06a24aca83d9591e1e45546bab652f448321854..535c416b47f7f4140dcba16f4e00e4a4f17af8d5 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -109,131 +109,6 @@ long long int runner_hist_bins[ runner_hist_N ];
 #endif
 
 
-/**
- * @brief Compute the 'interaction' between two particles.
- *
- * @param r2 the inter-particle radius squared.
- * @param hi the screening distance of the ith particle.
- * @param hj the screening distance of the jth particle.
- * @param io Pointer to where to store the interaction of the ith particle.
- * @param jo Pointer to where to store the interaction of the ith particle.
- */
- 
-__attribute__ ((always_inline)) INLINE void iact_nopart ( float r2 , float hi , float hj , float *force_i , float *force_j , int *count_i , int *count_j ) {
-
-    #define  KERNEL_COEFF_1  2.546479089470f
-    #define  KERNEL_COEFF_2  15.278874536822f
-    #define  KERNEL_COEFF_3  45.836623610466f
-    #define  KERNEL_COEFF_4  30.557749073644f
-    #define  KERNEL_COEFF_5  5.092958178941f
-    #define  KERNEL_COEFF_6  (-15.278874536822f)
-    #define  NORM_COEFF      4.188790204786f
-
-    float r = sqrtf( r2 );
-    float ui, uj, wi, wj;
-    
-    if ( r2 < hi*hi && !( force_i == NULL && count_i == NULL ) ) {
-        
-        ui = r / hi;
-        if ( ui < 0.5 )
-            wi = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (ui - 1.0f) * ui * ui;
-        else
-            wi = KERNEL_COEFF_5 * (1.0f - ui) * (1.0f - ui) * (1.0 - ui);
-        if ( force_i != NULL )
-            *force_i += NORM_COEFF * wi;
-        if ( count_i != NULL )
-            *count_i += 1;
-        
-        }
-
-    if ( r2 < hj*hj && !( force_j == NULL && count_j == NULL ) ) {
-        
-        uj = r / hj;
-        if ( uj < 0.5 )
-            wj = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (uj - 1.0f) * uj * uj;
-        else
-            wj = KERNEL_COEFF_5 * (1.0f - uj) * (1.0f - uj) * (1.0 - uj);
-        if ( force_j != NULL )
-            *force_j += NORM_COEFF * wj;
-        if ( count_j != NULL )
-            *count_j += 1;
-            
-        }
-        
-    #ifdef HIST
-    if ( hi > hj )
-        runner_hist_hit( hi / hj );
-    else
-        runner_hist_hit( hj / hi );
-    #endif
-    
-    }
-    
-
-
-__attribute__ ((always_inline)) INLINE void iact ( float r2 , float hi , float hj , struct part *pi , struct part *pj ) {
-
-    #define  KERNEL_COEFF_1  2.546479089470f
-    #define  KERNEL_COEFF_2  15.278874536822f
-    #define  KERNEL_COEFF_3  45.836623610466f
-    #define  KERNEL_COEFF_4  30.557749073644f
-    #define  KERNEL_COEFF_5  5.092958178941f
-    #define  KERNEL_COEFF_6  (-15.278874536822f)
-    #define  NORM_COEFF      4.188790204786f
-
-    float r = sqrtf( r2 );
-    float ui, uj, wi, wj;
-    float ui_dh, uj_dh, wi_dh, wj_dh;
-    
-    if ( r2 < hi*hi && pi != NULL ) {
-        
-        ui = r / hi;
-        ui_dh = -r / hi / hi;
-        if ( ui < 0.5f ) {
-            wi = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (ui - 1.0f) * ui * ui;
-            wi_dh = KERNEL_COEFF_2 * ui_dh * ui * ui
-                  + 2 * KERNEL_COEFF_2 * (ui - 1.0f) * ui_dh * ui;
-            }
-        else {
-            wi = KERNEL_COEFF_5 * (1.0f - ui) * (1.0f - ui) * (1.0f - ui);
-            wi_dh = -3 * KERNEL_COEFF_5 * ui_dh * (1.0f - ui) * (1.0f - ui);
-            }
-        pi->count += NORM_COEFF * wi;
-        pi->count_dh += NORM_COEFF * wi_dh;
-        pi->icount += 1;
-        
-        }
-
-    if ( r2 < hj*hj && pj != NULL ) {
-        
-        uj = r / hj;
-        uj_dh = -r / hj / hj;
-        if ( uj < 0.5f ) {
-            wj = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (uj - 1.0f) * uj * uj;
-            wj_dh = KERNEL_COEFF_2 * uj_dh * uj * uj
-                  + 2 * KERNEL_COEFF_2 * (uj - 1.0f) * uj_dh * uj;
-            }
-        else {
-            wj = KERNEL_COEFF_5 * (1.0f - uj) * (1.0f - uj) * (1.0f - uj);
-            wj_dh = -3 * KERNEL_COEFF_5 * uj_dh * (1.0f - uj) * (1.0f - uj);
-            }
-        pj->count += NORM_COEFF * wj;
-        pj->count_dh += NORM_COEFF * wj_dh;
-        pj->icount += 1;
-            
-        }
-        
-    #ifdef HIST
-    if ( hi > hj )
-        runner_hist_hit( hi / hj );
-    else
-        runner_hist_hit( hj / hi );
-    #endif
-    
-    }
-    
-
-
 /* A struct representing a runner's thread and its data. */
 struct runner_thread {
 
@@ -280,7 +155,8 @@ struct runner {
 
 /* Function prototypes. */
 void runner_run ( struct runner *r , int sort_queues );
-void runner_dopair ( struct runner_thread *rt , struct cell *ci , struct cell *cj );
-void runner_doself ( struct runner_thread *rt , struct cell *c );
+void runner_dopair_density ( struct runner_thread *rt , struct cell *ci , struct cell *cj );
+void runner_doself_density ( struct runner_thread *rt , struct cell *c );
+void runner_dosub_density ( struct runner_thread *rt , struct cell *ci , struct cell *cj , int flags );
 void runner_dosort ( struct runner_thread *rt , struct cell *c , int flag );
 void runner_init ( struct runner *r , struct space *s , int nr_threads , int nr_queues , int policy );
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
new file mode 100644
index 0000000000000000000000000000000000000000..a034e9d8641791b26fd2a444bf4204732faf8cd2
--- /dev/null
+++ b/src/runner_doiact.h
@@ -0,0 +1,543 @@
+/*******************************************************************************
+ * This file is part of GadgetSMP.
+ * Coypright (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/>.
+ * 
+ ******************************************************************************/
+
+
+/* Before including this file, define FUNCTION, which is the
+   name of the interaction function. This creates the interaction functions
+   runner_dopair_FUNCTION, runner_dopair_FUNCTION_naive, runner_doself_FUNCTION,
+   and runner_dosub_FUNCTION calling the pairwise interaction function
+   runner_iact_FUNCTION. */
+
+#define PASTE(x,y) x ## _ ## y
+
+#define DOPAIR2(f) PASTE(runner_dopair,f)
+#define DOPAIR DOPAIR2(FUNCTION)
+
+#define DOPAIR_NAIVE2(f) PASTE(runner_dopair_naive,f)
+#define DOPAIR_NAIVE DOPAIR_NAIVE2(FUNCTION)
+
+#define DOSELF2(f) PASTE(runner_doself,f)
+#define DOSELF DOSELF2(FUNCTION)
+
+#define DOSUB2(f) PASTE(runner_dosub,f)
+#define DOSUB DOSUB2(FUNCTION)
+
+#define IACT2(f) PASTE(runner_iact,f)
+#define IACT IACT2(FUNCTION)
+
+
+/**
+ * @brief Compute the interactions between a cell pair.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ */
+ 
+void DOPAIR_NAIVE ( struct runner_thread *rt , struct cell *ci , struct cell *cj ) {
+
+    struct runner *r = rt->r;
+    int pid, pjd, k, count_i = ci->count, count_j = cj->count;
+    double shift[3] = { 0.0 , 0.0 , 0.0 };
+    struct part *pi, *pj, *parts_i = ci->parts, *parts_j = cj->parts;
+    double dx[3], pix[3], hi, hi2, r2;
+    TIMER_TIC
+    
+    /* Get the relative distance between the pairs, wrapping. */
+    for ( k = 0 ; k < 3 ; k++ ) {
+        if ( cj->loc[k] - ci->loc[k] < -r->s->dim[k]/2 )
+            shift[k] = r->s->dim[k];
+        else if ( cj->loc[k] - ci->loc[k] > r->s->dim[k]/2 )
+            shift[k] = -r->s->dim[k];
+        }
+        
+    /* printf( "runner_dopair_naive: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts and shift = [ %g %g %g ].\n" ,
+        ci->loc[0] , ci->loc[1] , ci->loc[2] , cj->loc[0] , cj->loc[1] , cj->loc[2] ,
+        ci->count , cj->count , shift[0] , shift[1] , shift[2] ); fflush(stdout);
+    tic = getticks(); */
+    
+    /* Loop over the parts in ci. */
+    for ( pid = 0 ; pid < count_i ; pid++ ) {
+    
+        /* Get a hold of the ith part in ci. */
+        pi = &parts_i[ pid ];
+        for ( k = 0 ; k < 3 ; k++ )
+            pix[k] = pi->x[k] - shift[k];
+        hi = pi->h;
+        hi2 = hi * hi;
+        
+        /* Loop over the parts in cj. */
+        for ( pjd = 0 ; pjd < count_j ; pjd++ ) {
+        
+            /* Get a pointer to the jth particle. */
+            pj = &parts_j[ pjd ];
+        
+            /* Compute the pairwise distance. */
+            r2 = 0.0;
+            for ( k = 0 ; k < 3 ; k++ ) {
+                dx[k] = pix[k] - pj->x[k];
+                r2 += dx[k]*dx[k];
+                }
+                
+            /* Hit or miss? */
+            if ( r2 < hi2 || r2 < pj->h*pj->h ) {
+            
+                IACT( r2 , hi , pj->h , pi , pj );
+            
+                }
+        
+            } /* loop over the parts in cj. */
+    
+        } /* loop over the parts in ci. */
+        
+    #ifdef TIMER_VERBOSE
+        printf( "runner_dopair_naive[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f) took %.3f ms.\n" , rt->id , count_i , count_j , ci->depth , ci->r_max , cj->r_max , ((double)TIMER_TOC(runner_timer_dopair)) / CPU_TPS * 1000 );
+    #else
+        TIMER_TOC(runner_timer_dopair);
+    #endif
+
+
+    }
+
+
+/**
+ * @brief Compute the interactions between a cell pair.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ */
+ 
+void DOPAIR ( struct runner_thread *rt , struct cell *ci , struct cell *cj ) {
+
+    struct runner *r = rt->r;
+    int pid, pjd, k, sid;
+    double rshift, shift[3] = { 0.0 , 0.0 , 0.0 };
+    struct cell *temp;
+    struct entry *sort_i, *sort_j;
+    struct part *pi, *pj, *parts_i, *parts_j;
+    double dx[3], pix[3], pjx[3], hi, hi2, hj, hj2, r2, di, dj;
+    double hi_max, hj_max, di_max, dj_min;
+    int count_i, count_j;
+    TIMER_TIC
+    
+    /* Get the relative distance between the pairs, wrapping. */
+    for ( k = 0 ; k < 3 ; k++ ) {
+        if ( cj->loc[k] - ci->loc[k] < -r->s->dim[k]/2 )
+            shift[k] = r->s->dim[k];
+        else if ( cj->loc[k] - ci->loc[k] > r->s->dim[k]/2 )
+            shift[k] = -r->s->dim[k];
+        }
+        
+    /* Get the sorting index. */
+    for ( sid = 0 , k = 0 ; k < 3 ; k++ )
+        sid = 3*sid + ( (cj->loc[k] - ci->loc[k] + shift[k] < 0) ? 0 : (cj->loc[k] - ci->loc[k] + shift[k] > 0) ? 2 : 1 );
+
+    /* Switch the cells around? */
+    if ( runner_flip[sid] ) {
+        temp = ci; ci = cj; cj = temp;
+        for ( k = 0 ; k < 3 ; k++ )
+            shift[k] = -shift[k];
+        }
+    sid = sortlistID[sid];
+    
+    /* Get the cutoff shift. */
+    for ( rshift = 0.0 , k = 0 ; k < 3 ; k++ )
+        rshift += shift[k]*runner_shift[ 3*sid + k ];
+        
+    /* printf( "runner_dopair: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts and shift = [ %g %g %g ].\n" ,
+        ci->loc[0] , ci->loc[1] , ci->loc[2] , cj->loc[0] , cj->loc[1] , cj->loc[2] ,
+        ci->count , cj->count , shift[0] , shift[1] , shift[2] ); fflush(stdout); */
+    /* for ( hi = 0 , k = 0 ; k < ci->count ; k++ )
+        hi += ci->parts[k].r;
+    for ( hj = 0 , k = 0 ; k < cj->count ; k++ )
+        hj += cj->parts[k].r;
+    printf( "runner_dopair: avg. radii %g/%g for h=%g at depth=%i.\n" , hi/ci->count , hj/cj->count , ci->h[0] , ci->depth ); fflush(stdout); */
+    
+    /* Pick-out the sorted lists. */
+    sort_i = &ci->sort[ sid*(ci->count + 1) ];
+    sort_j = &cj->sort[ sid*(cj->count + 1) ];
+    
+    /* Get some other useful values. */
+    hi_max = ci->r_max - rshift; hj_max = cj->r_max - rshift;
+    count_i = ci->count; count_j = cj->count;
+    parts_i = ci->parts; parts_j = cj->parts;
+    di_max = sort_i[count_i-1].d - rshift;
+    dj_min = sort_j[0].d;
+    
+    /* if ( ci->split && cj->split && sid == 4 )
+        printf( "boing!\n" ); */
+    
+    /* Loop over the parts in ci. */
+    for ( pid = count_i-1 ; pid >= 0 && sort_i[pid].d + hi_max > dj_min ; pid-- ) {
+    
+        /* Get a hold of the ith part in ci. */
+        pi = &parts_i[ sort_i[ pid ].i ];
+        hi = pi->h;
+        di = sort_i[pid].d + hi - rshift;
+        if ( di < dj_min )
+            continue;
+            
+        hi2 = pi->h * pi->h;
+        for ( k = 0 ; k < 3 ; k++ )
+            pix[k] = pi->x[k] - shift[k];
+        
+        /* Loop over the parts in cj. */
+        for ( pjd = 0 ; pjd < count_j && sort_j[pjd].d < di ; pjd++ ) {
+        
+            /* Get a pointer to the jth particle. */
+            pj = &parts_j[ sort_j[pjd].i ];
+        
+            /* Compute the pairwise distance. */
+            r2 = 0.0;
+            for ( k = 0 ; k < 3 ; k++ ) {
+                dx[k] = pix[k] - pj->x[k];
+                r2 += dx[k]*dx[k];
+                }
+                
+            /* Hit or miss? */
+            if ( r2 < hi2 ) {
+            
+                IACT( r2 , hi , pj->h , pi , pj );
+            
+                }
+        
+            } /* loop over the parts in cj. */
+    
+        } /* loop over the parts in ci. */
+        
+    /* printf( "runner_dopair: first half took %.3f ms...\n" , ((double)(getticks() - tic)) / CPU_TPS * 1000 );
+    tic = getticks(); */
+
+    /* Loop over the parts in cj. */
+    for ( pjd = 0 ; pjd < count_j && sort_j[pjd].d - hj_max < di_max ; pjd++ ) {
+    
+        /* Get a hold of the jth part in cj. */
+        pj = &parts_j[ sort_j[ pjd ].i ];
+        hj = pj->h;
+        dj = sort_j[pjd].d - hj - rshift;
+        if ( dj > di_max )
+            continue;
+            
+        for ( k = 0 ; k < 3 ; k++ )
+            pjx[k] = pj->x[k] + shift[k];
+        hj2 = pj->h * pj->h;
+        
+        /* Loop over the parts in ci. */
+        for ( pid = count_i-1 ; pid >= 0 && sort_i[pid].d > dj ; pid-- ) {
+        
+            /* Get a pointer to the jth particle. */
+            pi = &parts_i[ sort_i[pid].i ];
+            
+            /* Compute the pairwise distance. */
+            r2 = 0.0;
+            for ( k = 0 ; k < 3 ; k++ ) {
+                dx[k] = pi->x[k] - pjx[k];
+                r2 += dx[k]*dx[k];
+                }
+                
+            /* Hit or miss? */
+            if ( r2 < hj2 && r2 > pi->h*pi->h ) {
+            
+                IACT( r2 , pi->h , hj , pi , pj );
+            
+                }
+        
+            } /* loop over the parts in cj. */
+    
+        } /* loop over the parts in ci. */
+
+    #ifdef TIMER_VERBOSE
+        printf( "runner_dopair[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f, h=%.3f) took %.3f ms.\n" , rt->id , count_i , count_j , ci->depth , ci->r_max , cj->r_max , fmax(ci->h[0],fmax(ci->h[1],ci->h[2])) , ((double)(TIMER_TOC(runner_timer_dopair))) / CPU_TPS * 1000 );
+    #else
+        TIMER_TOC(runner_timer_dopair);
+    #endif
+
+    }
+
+
+/**
+ * @brief Compute the cell self-interaction.
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ */
+
+void DOSELF ( struct runner_thread *rt , struct cell *c ) {
+
+    int k, pid, pjd, count = c->count;
+    double pix[3], dx[3], hi, hi2, r2;
+    struct part *pi, *pj, *parts = c->parts;
+    TIMER_TIC
+    
+    if ( c->split )
+        error( "Split cell should not have self-interactions." );
+    
+    /* Loop over the particles in the cell. */
+    for ( pid = 0 ; pid < count ; pid++ ) {
+    
+        /* Get a pointer to the ith particle. */
+        pi = &parts[pid];
+    
+        /* Get the particle position and radius. */
+        for ( k = 0 ; k < 3 ; k++ )
+            pix[k] = pi->x[k];
+        hi = pi->h;
+        hi2 = hi * hi;
+            
+        /* Loop over the other particles .*/
+        for ( pjd = pid+1 ; pjd < count ; pjd++ ) {
+        
+            /* Get a pointer to the jth particle. */
+            pj = &parts[pjd];
+        
+            /* Compute the pairwise distance. */
+            r2 = 0.0;
+            for ( k = 0 ; k < 3 ; k++ ) {
+                dx[k] = pix[k] - pj->x[k];
+                r2 += dx[k]*dx[k];
+                }
+                
+            /* Hit or miss? */
+            if ( r2 < hi2 || r2 < pj->h*pj->h ) {
+            
+                IACT( r2 , hi , pj->h , pi , pj );
+            
+                }
+        
+            } /* loop over all other particles. */
+    
+        } /* loop over all particles. */
+
+    #ifdef TIMER_VERBOSE
+        printf( "runner_doself[%02i]: %i parts at depth %i took %.3f ms.\n" , rt->id , count , c->depth , ((double)TIMER_TOC(runner_timer_doself)) / CPU_TPS * 1000 );
+    #else
+        TIMER_TOC(runner_timer_doself);
+    #endif
+
+    }
+
+
+/**
+ * @brief Compute grouped sub-cell interactions
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ */
+
+void DOSUB ( struct runner_thread *rt , struct cell *ci , struct cell *cj , int flags ) {
+
+    int j, k;
+
+    TIMER_TIC
+    
+    /* Different types of flags. */
+    switch ( flags ) {
+    
+        /* Regular sub-cell interactions of a single cell. */
+        case 0:
+            for ( j = 0 ; j < 7 ; j++ )
+                for ( k = j + 1 ; k < 8 ; k++ )
+                    if ( ci->progeny[j] != NULL && ci->progeny[k] != NULL )
+                        DOPAIR( rt , ci->progeny[j] , ci->progeny[k] );
+            break;
+            
+        case 1: /* (  1 ,  1 ,  0 ) */
+            if ( ci->progeny[6] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[0] );
+            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[1] );
+            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[0] );
+            if ( ci->progeny[7] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[1] );
+            break;
+    
+        case 3: /* (  1 ,  0 ,  1 ) */
+            if ( ci->progeny[5] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[0] );
+            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[2] );
+            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[0] );
+            if ( ci->progeny[7] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[2] );
+            break;
+                    
+        case 4: /* (  1 ,  0 ,  0 ) */
+            if ( ci->progeny[4] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[0] );
+            if ( ci->progeny[4] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[1] );
+            if ( ci->progeny[4] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[2] );
+            if ( ci->progeny[4] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[3] );
+            if ( ci->progeny[5] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[0] );
+            if ( ci->progeny[5] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[1] );
+            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[2] );
+            if ( ci->progeny[5] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[3] );
+            if ( ci->progeny[6] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[0] );
+            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[1] );
+            if ( ci->progeny[6] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[2] );
+            if ( ci->progeny[6] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[3] );
+            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[0] );
+            if ( ci->progeny[7] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[1] );
+            if ( ci->progeny[7] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[2] );
+            if ( ci->progeny[7] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[3] );
+            break;
+            
+        case 5: /* (  1 ,  0 , -1 ) */
+            if ( ci->progeny[4] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[1] );
+            if ( ci->progeny[4] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[3] );
+            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[1] );
+            if ( ci->progeny[6] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[3] );
+            break;
+                    
+        case 7: /* (  1 , -1 ,  0 ) */
+            if ( ci->progeny[4] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[2] );
+            if ( ci->progeny[4] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[4] , cj->progeny[3] );
+            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[2] );
+            if ( ci->progeny[5] != NULL && cj->progeny[3] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[3] );
+            break;
+                    
+        case 9: /* (  0 ,  1 ,  1 ) */
+            if ( ci->progeny[3] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[0] );
+            if ( ci->progeny[3] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[4] );
+            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[0] );
+            if ( ci->progeny[7] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[4] );
+            break;
+                    
+        case 10: /* (  0 ,  1 ,  0 ) */
+            if ( ci->progeny[2] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[2] , cj->progeny[0] );
+            if ( ci->progeny[2] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[2] , cj->progeny[1] );
+            if ( ci->progeny[2] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[2] , cj->progeny[4] );
+            if ( ci->progeny[2] != NULL && cj->progeny[5] != NULL )
+                DOPAIR( rt , ci->progeny[2] , cj->progeny[5] );
+            if ( ci->progeny[3] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[0] );
+            if ( ci->progeny[3] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[1] );
+            if ( ci->progeny[3] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[4] );
+            if ( ci->progeny[3] != NULL && cj->progeny[5] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[5] );
+            if ( ci->progeny[6] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[0] );
+            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[1] );
+            if ( ci->progeny[6] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[4] );
+            if ( ci->progeny[6] != NULL && cj->progeny[5] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[5] );
+            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[0] );
+            if ( ci->progeny[7] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[1] );
+            if ( ci->progeny[7] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[4] );
+            if ( ci->progeny[7] != NULL && cj->progeny[5] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[5] );
+            break;
+                    
+        case 11: /* (  0 ,  1 , -1 ) */
+            if ( ci->progeny[2] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[2] , cj->progeny[1] );
+            if ( ci->progeny[2] != NULL && cj->progeny[5] != NULL )
+                DOPAIR( rt , ci->progeny[2] , cj->progeny[5] );
+            if ( ci->progeny[6] != NULL && cj->progeny[1] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[1] );
+            if ( ci->progeny[6] != NULL && cj->progeny[5] != NULL )
+                DOPAIR( rt , ci->progeny[6] , cj->progeny[5] );
+            break;
+                    
+        case 12: /* (  0 ,  0 ,  1 ) */
+            if ( ci->progeny[1] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[1] , cj->progeny[0] );
+            if ( ci->progeny[1] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[1] , cj->progeny[2] );
+            if ( ci->progeny[1] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[1] , cj->progeny[4] );
+            if ( ci->progeny[1] != NULL && cj->progeny[6] != NULL )
+                DOPAIR( rt , ci->progeny[1] , cj->progeny[6] );
+            if ( ci->progeny[3] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[0] );
+            if ( ci->progeny[3] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[2] );
+            if ( ci->progeny[3] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[4] );
+            if ( ci->progeny[3] != NULL && cj->progeny[6] != NULL )
+                DOPAIR( rt , ci->progeny[3] , cj->progeny[6] );
+            if ( ci->progeny[5] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[0] );
+            if ( ci->progeny[5] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[2] );
+            if ( ci->progeny[5] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[4] );
+            if ( ci->progeny[5] != NULL && cj->progeny[6] != NULL )
+                DOPAIR( rt , ci->progeny[5] , cj->progeny[6] );
+            if ( ci->progeny[7] != NULL && cj->progeny[0] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[0] );
+            if ( ci->progeny[7] != NULL && cj->progeny[2] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[2] );
+            if ( ci->progeny[7] != NULL && cj->progeny[4] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[4] );
+            if ( ci->progeny[7] != NULL && cj->progeny[6] != NULL )
+                DOPAIR( rt , ci->progeny[7] , cj->progeny[6] );
+            break;
+                
+        }
+    
+
+    #ifdef TIMER_VERBOSE
+        printf( "runner_dosub[%02i]: flags=%i at depth %i took %.3f ms.\n" , rt->id , flags , ci->depth , ((double)TIMER_TOC(runner_timer_dosub)) / CPU_TPS * 1000 );
+    #else
+        TIMER_TOC(runner_timer_dosub);
+    #endif
+
+    }
+
+
diff --git a/src/runner_iact.h b/src/runner_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..d34296a5d1c79e9a0bea8f0cf3bdb14f320343c5
--- /dev/null
+++ b/src/runner_iact.h
@@ -0,0 +1,98 @@
+/*******************************************************************************
+ * This file is part of GadgetSMP.
+ * Coypright (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/>.
+ * 
+ ******************************************************************************/
+
+
+
+__attribute__ ((always_inline)) INLINE void runner_iact_density ( float r2 , float hi , float hj , struct part *pi , struct part *pj ) {
+
+    #define  KERNEL_COEFF_1  2.546479089470f
+    #define  KERNEL_COEFF_2  15.278874536822f
+    #define  KERNEL_COEFF_3  45.836623610466f
+    #define  KERNEL_COEFF_4  30.557749073644f
+    #define  KERNEL_COEFF_5  5.092958178941f
+    #define  KERNEL_COEFF_6  (-15.278874536822f)
+    #define  NORM_COEFF      4.188790204786f
+
+    float r = sqrtf( r2 );
+    float ui, uj, wi, wj;
+    float ui_dh, uj_dh, wi_dh, wj_dh;
+    
+    if ( r2 < hi*hi && pi != NULL ) {
+        
+        ui = r / hi;
+        ui_dh = -r / hi / hi;
+        if ( ui < 0.5f ) {
+            wi = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (ui - 1.0f) * ui * ui;
+            wi_dh = KERNEL_COEFF_2 * ui_dh * ui * ui
+                  + 2 * KERNEL_COEFF_2 * (ui - 1.0f) * ui_dh * ui;
+            }
+        else {
+            wi = KERNEL_COEFF_5 * (1.0f - ui) * (1.0f - ui) * (1.0f - ui);
+            wi_dh = -3 * KERNEL_COEFF_5 * ui_dh * (1.0f - ui) * (1.0f - ui);
+            }
+        pi->rho += NORM_COEFF * wi;
+        pi->rho_dh += NORM_COEFF * wi_dh;
+        pi->icount += 1;
+        
+        }
+
+    if ( r2 < hj*hj && pj != NULL ) {
+        
+        uj = r / hj;
+        uj_dh = -r / hj / hj;
+        if ( uj < 0.5f ) {
+            wj = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (uj - 1.0f) * uj * uj;
+            wj_dh = KERNEL_COEFF_2 * uj_dh * uj * uj
+                  + 2 * KERNEL_COEFF_2 * (uj - 1.0f) * uj_dh * uj;
+            }
+        else {
+            wj = KERNEL_COEFF_5 * (1.0f - uj) * (1.0f - uj) * (1.0f - uj);
+            wj_dh = -3 * KERNEL_COEFF_5 * uj_dh * (1.0f - uj) * (1.0f - uj);
+            }
+        pj->rho += NORM_COEFF * wj;
+        pj->rho_dh += NORM_COEFF * wj_dh;
+        pj->icount += 1;
+            
+        }
+        
+    #ifdef HIST
+    if ( hi > hj )
+        runner_hist_hit( hi / hj );
+    else
+        runner_hist_hit( hj / hi );
+    #endif
+    
+    }
+    
+
+
+__attribute__ ((always_inline)) INLINE void runner_iact_force ( float r2 , float hi , float hj , struct part *pi , struct part *pj ) {
+
+    #define  KERNEL_COEFF_1  2.546479089470f
+    #define  KERNEL_COEFF_2  15.278874536822f
+    #define  KERNEL_COEFF_3  45.836623610466f
+    #define  KERNEL_COEFF_4  30.557749073644f
+    #define  KERNEL_COEFF_5  5.092958178941f
+    #define  KERNEL_COEFF_6  (-15.278874536822f)
+    #define  NORM_COEFF      4.188790204786f
+    
+    }
+    
+
+
diff --git a/src/space.c b/src/space.c
index c8795a0b9961f231fe9a2a96468f3a77d2222a95..54037106a8880dd8f6d7916fd12711df634f6628 100644
--- a/src/space.c
+++ b/src/space.c
@@ -80,7 +80,7 @@ const int sortlistID[27] = {
     
     
 /**
- * @brief Mapping function to draw a specific cell (gnuplot).
+ * @brief Mapping function to free the sorted indices buffers.
  */
 
 void space_map_clearsort ( struct cell *c , void *data ) {
@@ -93,6 +93,20 @@ void space_map_clearsort ( struct cell *c , void *data ) {
     }
 
 
+/**
+ * @brief Mapping function to append a ghost task to each cell.
+ */
+
+void space_map_mkghosts ( struct cell *c , void *data ) {
+
+    struct space *s = (struct space *)data;
+
+    /* Make the ghost task. */
+    c->ghost = space_addtask( s , task_type_ghost , task_subtype_none , 0 , 0 , c , NULL , NULL , 0 , NULL , 0 );
+
+    }
+
+
 /**
  * @brief Get a task free of dependencies and conflicts.
  *
@@ -124,15 +138,15 @@ struct task *space_gettask ( struct space *s ) {
             
             /* Different criteria for different types. */
             switch ( res->type ) {
-                case tid_self:
+                case task_type_self:
                     if ( res->ci->lock || res->ci->wait > 0 )
                         continue;
                     break;
-                case tid_pair:
+                case task_type_pair:
                     if ( res->ci->lock || res->cj->lock || res->ci->wait || res->cj->wait )
                         continue;
                     break;
-                case tid_sort:
+                case task_type_sort:
                     if ( res->ci->lock )
                         continue;
                     break;
@@ -160,7 +174,7 @@ struct task *space_gettask ( struct space *s ) {
             s->tasks_ind[ s->next_task ] = tid;
             
             /* Lock the cells, if needed. */
-            if ( s->tasks[tid].type != tid_sort ) {
+            if ( s->tasks[tid].type != task_type_sort ) {
                 for ( c = res->ci ; c != NULL ; c = c->parent )
                     __sync_fetch_and_add( &c->lock , 1 );
                 for ( c = res->cj ; c != NULL ; c = c->parent )
@@ -226,9 +240,10 @@ void space_map_parts ( struct space *s , void (*fun)( struct part *p , struct ce
  * @brief Map a function to all particles in a aspace.
  *
  * @param s The #space we are working in.
+ * @param full Map to all cells, including cells with sub-cells.
  */
  
-void space_map_cells ( struct space *s , void (*fun)( struct cell *c , void *data ) , void *data ) {
+void space_map_cells ( struct space *s , int full , void (*fun)( struct cell *c , void *data ) , void *data ) {
 
     int i;
 
@@ -237,11 +252,11 @@ void space_map_cells ( struct space *s , void (*fun)( struct cell *c , void *dat
         int k;
         
         /* No progeny? */
-        if ( !c->split )
+        if ( full || !c->split )
             fun( c , data );
                 
-        /* Otherwise, recurse. */
-        else
+        /* Recurse. */
+        if ( c->split )
             for ( k = 0 ; k < 8 ; k++ )
                 if ( c->progeny[k] != NULL )
                     rec_map( c->progeny[k] );
@@ -261,12 +276,13 @@ void space_map_cells ( struct space *s , void (*fun)( struct cell *c , void *dat
  * @param s The #space we are working in.
  */
  
-struct task *space_addtask ( struct space *s , int type , int flags , int wait , struct cell *ci , struct cell *cj , struct task *unlock_tasks[] , int nr_unlock_tasks , struct cell *unlock_cells[] , int nr_unlock_cells ) {
+struct task *space_addtask ( struct space *s , int type , int subtype , int flags , int wait , struct cell *ci , struct cell *cj , struct task *unlock_tasks[] , int nr_unlock_tasks , struct cell *unlock_cells[] , int nr_unlock_cells ) {
 
     struct task *t = &s->tasks[ s->nr_tasks ];
     
     /* Copy the data. */
     t->type = type;
+    t->subtype = subtype;
     t->flags = flags;
     t->wait = wait;
     t->ci = ci;
@@ -308,7 +324,7 @@ void space_splittasks ( struct space *s ) {
         t = &s->tasks[tid];
     
         /* If this task isn't a pair, i'm not interested. */
-        if ( t->type != tid_pair )
+        if ( t->type != task_type_pair )
             continue;
             
         /* Get a handle on the cells involved. */
@@ -364,7 +380,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[6]->split && !ci->progeny[7]->split &&
                          !cj->progeny[0]->split && !cj->progeny[1]->split ) {
-                        t->type = tid_sub; t->flags = 1;
+                        t->type = task_type_sub; t->flags = 1;
                         task_addunlock( ci->progeny[6]->sorts[1] , t ); task_addunlock( cj->progeny[0]->sorts[1] , t );
                         task_addunlock( ci->progeny[7]->sorts[1] , t ); task_addunlock( cj->progeny[1]->sorts[1] , t );
                         task_addunlock( ci->progeny[6]->sorts[0] , t ); task_addunlock( cj->progeny[1]->sorts[0] , t );
@@ -373,11 +389,11 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[6]; t->cj = cj->progeny[0];
                         task_addunlock( ci->progeny[6]->sorts[1] , t ); task_addunlock( cj->progeny[0]->sorts[1] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[1] , t ); task_addunlock( cj->progeny[1]->sorts[1] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[0] , t ); task_addunlock( cj->progeny[1]->sorts[0] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[2] , t ); task_addunlock( cj->progeny[0]->sorts[2] , t );
                         }
                     ci->progeny[6]->nr_pairs += 2;
@@ -397,7 +413,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[5]->split && !ci->progeny[7]->split &&
                          !cj->progeny[0]->split && !cj->progeny[2]->split ) {
-                        t->type = tid_sub; t->flags = 3;
+                        t->type = task_type_sub; t->flags = 3;
                         task_addunlock( ci->progeny[5]->sorts[3] , t ); task_addunlock( cj->progeny[0]->sorts[3] , t );
                         task_addunlock( ci->progeny[7]->sorts[3] , t ); task_addunlock( cj->progeny[2]->sorts[3] , t );
                         task_addunlock( ci->progeny[5]->sorts[0] , t ); task_addunlock( cj->progeny[2]->sorts[0] , t );
@@ -406,11 +422,11 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[5]; t->cj = cj->progeny[0];
                         task_addunlock( ci->progeny[5]->sorts[3] , t ); task_addunlock( cj->progeny[0]->sorts[3] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[3] , t ); task_addunlock( cj->progeny[2]->sorts[3] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[0] , t ); task_addunlock( cj->progeny[2]->sorts[0] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[6] , t ); task_addunlock( cj->progeny[0]->sorts[6] , t );
                         }
                     ci->progeny[5]->nr_pairs += 2;
@@ -423,7 +439,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[4]->split && !ci->progeny[5]->split && !ci->progeny[6]->split && !ci->progeny[7]->split &&
                          !cj->progeny[0]->split && !cj->progeny[1]->split && !cj->progeny[2]->split && !cj->progeny[3]->split ) {
-                        t->type = tid_sub; t->flags = 4;
+                        t->type = task_type_sub; t->flags = 4;
                         task_addunlock( ci->progeny[4]->sorts[4] , t ); task_addunlock( cj->progeny[0]->sorts[4] , t );
                         task_addunlock( ci->progeny[5]->sorts[5] , t ); task_addunlock( cj->progeny[0]->sorts[5] , t );
                         task_addunlock( ci->progeny[6]->sorts[7] , t ); task_addunlock( cj->progeny[0]->sorts[7] , t );
@@ -444,35 +460,35 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[4]; t->cj = cj->progeny[0];
                         task_addunlock( ci->progeny[4]->sorts[4] , t ); task_addunlock( cj->progeny[0]->sorts[4] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[5] , t ); task_addunlock( cj->progeny[0]->sorts[5] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[7] , t ); task_addunlock( cj->progeny[0]->sorts[7] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[8] , t ); task_addunlock( cj->progeny[0]->sorts[8] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[4] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[4] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[4]->sorts[3] , t ); task_addunlock( cj->progeny[1]->sorts[3] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[4] , t ); task_addunlock( cj->progeny[1]->sorts[4] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[6] , t ); task_addunlock( cj->progeny[1]->sorts[6] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[7] , t ); task_addunlock( cj->progeny[1]->sorts[7] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[4] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[4] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[4]->sorts[1] , t ); task_addunlock( cj->progeny[2]->sorts[1] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[2] , t ); task_addunlock( cj->progeny[2]->sorts[2] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[4] , t ); task_addunlock( cj->progeny[2]->sorts[4] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[5] , t ); task_addunlock( cj->progeny[2]->sorts[5] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[4] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[4] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[4]->sorts[0] , t ); task_addunlock( cj->progeny[3]->sorts[0] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[1] , t ); task_addunlock( cj->progeny[3]->sorts[1] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[3] , t ); task_addunlock( cj->progeny[3]->sorts[3] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[4] , t ); task_addunlock( cj->progeny[3]->sorts[4] , t );
                         }
                     ci->progeny[4]->nr_pairs += 4;
@@ -489,7 +505,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[4]->split && !ci->progeny[6]->split &&
                          !cj->progeny[1]->split && !cj->progeny[3]->split ) {
-                        t->type = tid_sub; t->flags = 5;
+                        t->type = task_type_sub; t->flags = 5;
                         task_addunlock( ci->progeny[4]->sorts[5] , t ); task_addunlock( cj->progeny[1]->sorts[5] , t );
                         task_addunlock( ci->progeny[6]->sorts[5] , t ); task_addunlock( cj->progeny[3]->sorts[5] , t );
                         task_addunlock( ci->progeny[4]->sorts[2] , t ); task_addunlock( cj->progeny[3]->sorts[2] , t );
@@ -498,11 +514,11 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[4]; t->cj = cj->progeny[1];
                         task_addunlock( ci->progeny[4]->sorts[5] , t ); task_addunlock( cj->progeny[1]->sorts[5] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[5] , t ); task_addunlock( cj->progeny[3]->sorts[5] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[4] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[4] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[4]->sorts[2] , t ); task_addunlock( cj->progeny[3]->sorts[2] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[8] , t ); task_addunlock( cj->progeny[1]->sorts[8] , t );
                         }
                     ci->progeny[4]->nr_pairs += 2;
@@ -522,7 +538,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[4]->split && !ci->progeny[5]->split &&
                          !cj->progeny[2]->split && !cj->progeny[3]->split ) {
-                        t->type = tid_sub; t->flags = 7;
+                        t->type = task_type_sub; t->flags = 7;
                         task_addunlock( ci->progeny[4]->sorts[6] , t ); task_addunlock( cj->progeny[3]->sorts[6] , t );
                         task_addunlock( ci->progeny[5]->sorts[8] , t ); task_addunlock( cj->progeny[2]->sorts[8] , t );
                         task_addunlock( ci->progeny[4]->sorts[7] , t ); task_addunlock( cj->progeny[2]->sorts[7] , t );
@@ -531,11 +547,11 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[4]; t->cj = cj->progeny[3];
                         task_addunlock( ci->progeny[4]->sorts[6] , t ); task_addunlock( cj->progeny[3]->sorts[6] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[8] , t ); task_addunlock( cj->progeny[2]->sorts[8] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[4] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[4] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[4]->sorts[7] , t ); task_addunlock( cj->progeny[2]->sorts[7] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[3] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[3] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[7] , t ); task_addunlock( cj->progeny[3]->sorts[7] , t );
                         }
                     ci->progeny[4]->nr_pairs += 2;
@@ -555,7 +571,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[3]->split && !ci->progeny[7]->split &&
                          !cj->progeny[0]->split && !cj->progeny[4]->split ) {
-                        t->type = tid_sub; t->flags = 9;
+                        t->type = task_type_sub; t->flags = 9;
                         task_addunlock( ci->progeny[3]->sorts[9] , t ); task_addunlock( cj->progeny[0]->sorts[9] , t );
                         task_addunlock( ci->progeny[7]->sorts[9] , t ); task_addunlock( cj->progeny[4]->sorts[9] , t );
                         task_addunlock( ci->progeny[3]->sorts[0] , t ); task_addunlock( cj->progeny[4]->sorts[0] , t );
@@ -564,11 +580,11 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[3]; t->cj = cj->progeny[0];
                         task_addunlock( ci->progeny[3]->sorts[9] , t ); task_addunlock( cj->progeny[0]->sorts[9] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[9] , t ); task_addunlock( cj->progeny[4]->sorts[9] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[0] , t ); task_addunlock( cj->progeny[4]->sorts[0] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[8] , t ); task_addunlock( cj->progeny[0]->sorts[8] , t );
                         }
                     ci->progeny[3]->nr_pairs += 2;
@@ -580,7 +596,7 @@ void space_splittasks ( struct space *s ) {
                 case 10: /* (  0 ,  1 ,  0 ) */
                     if ( !ci->progeny[2]->split && !ci->progeny[3]->split && !ci->progeny[6]->split && !ci->progeny[7]->split &&
                          !cj->progeny[0]->split && !cj->progeny[1]->split && !cj->progeny[4]->split && !cj->progeny[5]->split ) {
-                        t->type = tid_sub; t->flags = 10;
+                        t->type = task_type_sub; t->flags = 10;
                         task_addunlock( ci->progeny[2]->sorts[10] , t ); task_addunlock( cj->progeny[0]->sorts[10] , t );
                         task_addunlock( ci->progeny[3]->sorts[11] , t ); task_addunlock( cj->progeny[0]->sorts[11] , t );
                         task_addunlock( ci->progeny[6]->sorts[7] , t ); task_addunlock( cj->progeny[0]->sorts[7] , t );
@@ -601,35 +617,35 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[2]; t->cj = cj->progeny[0];
                         task_addunlock( ci->progeny[2]->sorts[10] , t ); task_addunlock( cj->progeny[0]->sorts[10] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[11] , t ); task_addunlock( cj->progeny[0]->sorts[11] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[7] , t ); task_addunlock( cj->progeny[0]->sorts[7] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[6] , t ); task_addunlock( cj->progeny[0]->sorts[6] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[2] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[2] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[2]->sorts[9] , t ); task_addunlock( cj->progeny[1]->sorts[9] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[10] , t ); task_addunlock( cj->progeny[1]->sorts[10] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[8] , t ); task_addunlock( cj->progeny[1]->sorts[8] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[7] , t ); task_addunlock( cj->progeny[1]->sorts[7] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[2] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[2] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[2]->sorts[1] , t ); task_addunlock( cj->progeny[4]->sorts[1] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[2] , t ); task_addunlock( cj->progeny[4]->sorts[2] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[10] , t ); task_addunlock( cj->progeny[4]->sorts[10] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[11] , t ); task_addunlock( cj->progeny[4]->sorts[11] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[2] , cj->progeny[5] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[2] , cj->progeny[5] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[2]->sorts[0] , t ); task_addunlock( cj->progeny[5]->sorts[0] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[5] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[5] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[1] , t ); task_addunlock( cj->progeny[5]->sorts[1] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[5] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[5] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[9] , t ); task_addunlock( cj->progeny[5]->sorts[9] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[5] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[5] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[10] , t ); task_addunlock( cj->progeny[5]->sorts[10] , t );
                         }
                     ci->progeny[2]->nr_pairs += 4;
@@ -646,7 +662,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[2]->split && !ci->progeny[6]->split &&
                          !cj->progeny[1]->split && !cj->progeny[5]->split ) {
-                        t->type = tid_sub; t->flags = 11;
+                        t->type = task_type_sub; t->flags = 11;
                         task_addunlock( ci->progeny[2]->sorts[11] , t ); task_addunlock( cj->progeny[1]->sorts[11] , t );
                         task_addunlock( ci->progeny[6]->sorts[11] , t ); task_addunlock( cj->progeny[5]->sorts[11] , t );
                         task_addunlock( ci->progeny[2]->sorts[2] , t ); task_addunlock( cj->progeny[5]->sorts[2] , t );
@@ -655,11 +671,11 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[2]; t->cj = cj->progeny[1];
                         task_addunlock( ci->progeny[2]->sorts[11] , t ); task_addunlock( cj->progeny[1]->sorts[11] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[5] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[5] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[11] , t ); task_addunlock( cj->progeny[5]->sorts[11] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[2] , cj->progeny[5] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[2] , cj->progeny[5] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[2]->sorts[2] , t ); task_addunlock( cj->progeny[5]->sorts[2] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[6] , cj->progeny[1] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[6]->sorts[6] , t ); task_addunlock( cj->progeny[1]->sorts[6] , t );
                         }
                     ci->progeny[2]->nr_pairs += 2;
@@ -672,7 +688,7 @@ void space_splittasks ( struct space *s ) {
                     if ( space_dosub &&
                          !ci->progeny[1]->split && !ci->progeny[3]->split && !ci->progeny[5]->split && !ci->progeny[7]->split &&
                          !cj->progeny[0]->split && !cj->progeny[2]->split && !cj->progeny[4]->split && !cj->progeny[6]->split ) {
-                        t->type = tid_sub; t->flags = 12;
+                        t->type = task_type_sub; t->flags = 12;
                         task_addunlock( ci->progeny[1]->sorts[12] , t ); task_addunlock( cj->progeny[0]->sorts[12] , t );
                         task_addunlock( ci->progeny[3]->sorts[11] , t ); task_addunlock( cj->progeny[0]->sorts[11] , t );
                         task_addunlock( ci->progeny[5]->sorts[5] , t ); task_addunlock( cj->progeny[0]->sorts[5] , t );
@@ -693,35 +709,35 @@ void space_splittasks ( struct space *s ) {
                     else {
                         t->ci = ci->progeny[1]; t->cj = cj->progeny[0];
                         task_addunlock( ci->progeny[1]->sorts[12] , t ); task_addunlock( cj->progeny[0]->sorts[12] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[11] , t ); task_addunlock( cj->progeny[0]->sorts[11] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[5] , t ); task_addunlock( cj->progeny[0]->sorts[5] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[0] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[2] , t ); task_addunlock( cj->progeny[0]->sorts[2] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[1] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[1] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[1]->sorts[9] , t ); task_addunlock( cj->progeny[2]->sorts[9] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[12] , t ); task_addunlock( cj->progeny[2]->sorts[12] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[8] , t ); task_addunlock( cj->progeny[2]->sorts[8] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[2] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[2] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[5] , t ); task_addunlock( cj->progeny[2]->sorts[5] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[1] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[1] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[1]->sorts[3] , t ); task_addunlock( cj->progeny[4]->sorts[3] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[6] , t ); task_addunlock( cj->progeny[4]->sorts[6] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[12] , t ); task_addunlock( cj->progeny[4]->sorts[12] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[4] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[4] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[11] , t ); task_addunlock( cj->progeny[4]->sorts[11] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[1] , cj->progeny[6] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[1] , cj->progeny[6] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[1]->sorts[0] , t ); task_addunlock( cj->progeny[6]->sorts[0] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[3] , cj->progeny[6] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[3] , cj->progeny[6] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[3]->sorts[3] , t ); task_addunlock( cj->progeny[6]->sorts[3] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[5] , cj->progeny[6] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[5] , cj->progeny[6] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[5]->sorts[9] , t ); task_addunlock( cj->progeny[6]->sorts[9] , t );
-                        t = space_addtask( s , tid_pair , 0 , 0 , ci->progeny[7] , cj->progeny[6] , NULL , 0 , NULL , 0 );
+                        t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , ci->progeny[7] , cj->progeny[6] , NULL , 0 , NULL , 0 );
                         task_addunlock( ci->progeny[7]->sorts[12] , t ); task_addunlock( cj->progeny[6]->sorts[12] , t );
                         }
                     ci->progeny[1]->nr_pairs += 4;
@@ -758,7 +774,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
     int i, j, k, ii, jj, kk, iii, jjj, kkk, cid, cjd;
     int *cdim = s->cdim;
     int nr_tasks_old = s->nr_tasks;
-    struct task *t;
+    struct task *t, *t2;
     int pts[7][8] = { { -1 , 12 , 10 , 9 , 4 , 3 , 1 , 0 } ,
                       { -1 , -1 , 11 , 10 , 5 , 4 , 2 , 1 } ,
                       { -1 , -1 , -1 , 12 , 7 , 6 , 4 , 3 } , 
@@ -766,7 +782,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
                       { -1 , -1 , -1 , -1 , -1 , 12 , 10 , 9 } ,
                       { -1 , -1 , -1 , -1 , -1 , -1 , 11 , 10 } ,
                       { -1 , -1 , -1 , -1 , -1 , -1 , -1 , 12 } };
-    int counts[tid_count];
+    int counts[task_type_count];
 
     /* Recursive function to generate tasks in the cell tree. */
     void maketasks_rec ( struct cell *c , struct task *sort_up[] , int nr_sort_up , struct cell *parent ) {
@@ -783,15 +799,15 @@ void space_maketasks ( struct space *s , int do_sort ) {
         
             if ( do_sort ) {
                 if ( c->count < 1000 ) {
-                    sort[0] = space_addtask( s , tid_sort , 0x1fff , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    sort[0] = space_addtask( s , task_type_sort , task_subtype_none , 0x1fff , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
                     for ( k = 0 ; k < 13 ; k++ )
                         c->sorts[k] = sort[0];
                     nr_sort = 1;
                     }
                 else if ( c->count < 5000 ) {
-                    sort[0] = space_addtask( s , tid_sort , 0xf , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    sort[1] = space_addtask( s , tid_sort , 0xf0 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    sort[2] = space_addtask( s , tid_sort , 0x1f00 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    sort[0] = space_addtask( s , task_type_sort , task_subtype_none , 0xf , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    sort[1] = space_addtask( s , task_type_sort , task_subtype_none , 0xf0 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    sort[2] = space_addtask( s , task_type_sort , task_subtype_none , 0x1f00 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
                     for ( k = 0 ; k < 4 ; k++ )
                         c->sorts[k] = sort[0];
                     for ( k = 4 ; k < 8 ; k++ )
@@ -801,26 +817,26 @@ void space_maketasks ( struct space *s , int do_sort ) {
                     nr_sort = 3;
                     }
                 else {
-                    c->sorts[0] = sort[0] = space_addtask( s , tid_sort , 0x1 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[1] = sort[1] = space_addtask( s , tid_sort , 0x2 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[2] = sort[2] = space_addtask( s , tid_sort , 0x4 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[3] = sort[3] = space_addtask( s , tid_sort , 0x8 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[4] = sort[4] = space_addtask( s , tid_sort , 0x10 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[5] = sort[5] = space_addtask( s , tid_sort , 0x20 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[6] = sort[6] = space_addtask( s , tid_sort , 0x40 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[7] = sort[7] = space_addtask( s , tid_sort , 0x80 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[8] = sort[8] = space_addtask( s , tid_sort , 0x100 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[9] = sort[9] = space_addtask( s , tid_sort , 0x200 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[10] = sort[10] = space_addtask( s , tid_sort , 0x400 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[11] = sort[11] = space_addtask( s , tid_sort , 0x800 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
-                    c->sorts[12] = sort[12] = space_addtask( s , tid_sort , 0x1000 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[0] = sort[0] = space_addtask( s , task_type_sort , task_subtype_none , 0x1 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[1] = sort[1] = space_addtask( s , task_type_sort , task_subtype_none , 0x2 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[2] = sort[2] = space_addtask( s , task_type_sort , task_subtype_none , 0x4 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[3] = sort[3] = space_addtask( s , task_type_sort , task_subtype_none , 0x8 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[4] = sort[4] = space_addtask( s , task_type_sort , task_subtype_none , 0x10 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[5] = sort[5] = space_addtask( s , task_type_sort , task_subtype_none , 0x20 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[6] = sort[6] = space_addtask( s , task_type_sort , task_subtype_none , 0x40 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[7] = sort[7] = space_addtask( s , task_type_sort , task_subtype_none , 0x80 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[8] = sort[8] = space_addtask( s , task_type_sort , task_subtype_none , 0x100 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[9] = sort[9] = space_addtask( s , task_type_sort , task_subtype_none , 0x200 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[10] = sort[10] = space_addtask( s , task_type_sort , task_subtype_none , 0x400 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[11] = sort[11] = space_addtask( s , task_type_sort , task_subtype_none , 0x800 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
+                    c->sorts[12] = sort[12] = space_addtask( s , task_type_sort , task_subtype_none , 0x1000 , 0 , c , NULL , sort_up , nr_sort_up , NULL , 0 );
                     nr_sort = 13;
                     }
                 }
 
             /* Generate a self-interaction if not split. */
             if ( !c->split && c->count > 1 )
-                space_addtask( s , tid_self , 0 , 0 , c , NULL , NULL , 0 , NULL , 0 );
+                space_addtask( s , task_type_self , task_subtype_density , 0 , 0 , c , NULL , NULL , 0 , NULL , 0 );
                 
             }
             
@@ -840,7 +856,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
                     if ( c->progeny[j] != NULL && c->progeny[j]->count > 0 )
                         for ( k = j + 1 ; k < 8 ; k++ )
                             if ( c->progeny[k] != NULL && c->progeny[k]->count > 0 ) {
-                                t = space_addtask( s , tid_pair , 0 , 0 , c->progeny[j] , c->progeny[k] , NULL , 0 , NULL , 0 );
+                                t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , c->progeny[j] , c->progeny[k] , NULL , 0 , NULL , 0 );
                                 task_addunlock( c->progeny[j]->sorts[ pts[j][k] ] , t );
                                 task_addunlock( c->progeny[k]->sorts[ pts[j][k] ] , t );
                                 c->progeny[k]->nr_pairs += 1;
@@ -853,7 +869,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
             else {
             
                 /* Add the task. */
-                t = space_addtask( s , tid_sub , 0 , 0 , c , NULL , NULL , 0 , NULL , 0 );
+                t = space_addtask( s , task_type_sub , task_subtype_density , 0 , 0 , c , NULL , NULL , 0 , NULL , 0 );
                 
                 /* Make it depend on all the sorts of its progeny. */
                 for ( k = 0 ; k < 8 ; k++ )
@@ -904,7 +920,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
                                 continue;
                             if ( cid >= cjd )
                                 continue;
-                            t = space_addtask( s , tid_pair , 0 , 0 , &s->cells[cid] , &s->cells[cjd] , NULL , 0 , NULL , 0 );
+                            t = space_addtask( s , task_type_pair , task_subtype_density , 0 , 0 , &s->cells[cid] , &s->cells[cjd] , NULL , 0 , NULL , 0 );
                             task_addunlock( s->cells[cid].sorts[ sortlistID[ (kk+1) + 3*( (jj+1) + 3*(ii+1) ) ] ] , t );
                             task_addunlock( s->cells[cjd].sorts[ sortlistID[ (kk+1) + 3*( (jj+1) + 3*(ii+1) ) ] ] , t );
                             s->cells[cid].nr_pairs += 1;
@@ -916,6 +932,35 @@ void space_maketasks ( struct space *s , int do_sort ) {
 
     /* Split the tasks. */
     space_splittasks( s );
+    
+    /* Append a ghost task to each cell. */
+    space_map_cells( s , 1 , &space_map_mkghosts , s );
+    
+    /* Run through the tasks and make iacts for each density task. */
+    for ( k = 0 ; k < s->nr_tasks ; k++ ) {
+    
+        /* Get a pointer to the task. */
+        t = &s->tasks[k];
+        
+        /* Self-interaction? */
+        if ( t->type == task_type_self && t->subtype == task_subtype_density ) {
+            task_addunlock( t , t->ci->ghost );
+            t2 = space_addtask( s , task_type_self , task_subtype_force , 0 , 0 , t->ci , NULL , NULL , 0 , NULL , 0 );
+            task_addunlock( t->ci->ghost , t2 );
+            }
+            
+        /* Otherwise, pair interaction? */
+        else if ( t->type == task_type_pair && t->subtype == task_subtype_density ) {
+            if ( t->ci->ghost == NULL )
+                printf( "space_maketasks: cell at %lx has no ghost!\n" , (unsigned long int)t->ci );
+            task_addunlock( t , t->ci->ghost );
+            task_addunlock( t , t->cj->ghost );
+            t2 = space_addtask( s , task_type_pair , task_subtype_force , 0 , 0 , t->ci , t->cj , NULL , 0 , NULL , 0 );
+            task_addunlock( t->ci->ghost , t2 );
+            task_addunlock( t->cj->ghost , t2 );
+            }
+
+        }
         
     /* Did we already create indices? */
     if ( s->tasks_ind == NULL )
@@ -927,11 +972,11 @@ void space_maketasks ( struct space *s , int do_sort ) {
         for ( k = 0 ; k < s->nr_tasks ; k++ )
             s->tasks_ind[k] = k;
             
-    /* Remove sort tasks with no dependencies. */
+    /* Remove sort and ghost tasks with no dependencies. */
     for ( k = 0 ; k < s->nr_tasks ; k++ ) {
         t = &s->tasks[k];
-        if ( t->type == tid_sort && t->nr_unlock_tasks == 0 ) {
-            t->type = tid_none;
+        if ( ( t->type == task_type_sort || t->type == task_type_ghost ) && t->nr_unlock_tasks == 0 ) {
+            t->type = task_type_none;
             if ( t->ci->split )
                 for ( j = 0 ; j < 8 ; j++ )
                     if ( t->ci->progeny[j] != NULL && t->flags & ( 1 << j ) )
@@ -940,12 +985,12 @@ void space_maketasks ( struct space *s , int do_sort ) {
         }
             
     /* Count the number of each task type. */
-    for ( k = 0 ; k < tid_count ; k++ )
+    for ( k = 0 ; k < task_type_count ; k++ )
         counts[k] = 0;
     for ( k = 0 ; k < s->nr_tasks ; k++ )
         counts[ s->tasks[k].type ] += 1;
     printf( "space_maketasks: task counts are [ %s=%i" , taskID_names[0] , counts[0] );
-    for ( k = 1 ; k < tid_count ; k++ )
+    for ( k = 1 ; k < task_type_count ; k++ )
         printf( " %s=%i" , taskID_names[k] , counts[k] );
     printf( " ]\n" );
         
@@ -966,7 +1011,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
 void space_split ( struct space *s , struct cell *c ) {
 
     int k, count;
-    double r, r_limit, r_max = 0.0;
+    double h, h_limit, r_max = 0.0;
     struct cell *temp;
     
     /* Check the depth. */
@@ -974,15 +1019,15 @@ void space_split ( struct space *s , struct cell *c ) {
         s->maxdepth = c->depth;
     
     /* Set the minimum cutoff. */
-    r_limit = fmin( c->h[0] , fmin( c->h[1] , c->h[2] ) ) / 2;
+    h_limit = fmin( c->h[0] , fmin( c->h[1] , c->h[2] ) ) / 2;
     
     /* Count the particles below that. */
     for ( count = 0 , k = 0 ; k < c->count ; k++ ) {
-        r = c->parts[k].r;
-        if ( r <= r_limit )
+        h = c->parts[k].h;
+        if ( h <= h_limit )
             count += 1;
-        if ( r > r_max )
-            r_max = r;
+        if ( h > r_max )
+            r_max = h;
         }
     c->r_max = r_max;
             
@@ -1113,28 +1158,28 @@ struct cell *space_getcell ( struct space *s ) {
  */
 
 
-void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_max ) {
+void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_cells ) {
 
     int i, j, k;
     int nr_cells, cdim[3];
-    double r_min, r_max, h[3], ih[3];
+    double h_min, h_max, h[3], ih[3];
     struct cell *c, *cells;
     struct part *parts_new, *finger;
     
     
     /* Get the minimum and maximum cutoff radii. */
-    r_min = parts[0].r; r_max = r_min;
+    h_min = parts[0].h; h_max = h_min;
     for ( k = 1 ; k < N ; k++ )
-        if ( parts[k].r < r_min )
-            r_min = parts[k].r;
-        else if ( parts[k].r > r_max )
-            r_max = parts[k].r;
+        if ( parts[k].h < h_min )
+            h_min = parts[k].h;
+        else if ( parts[k].h > h_max )
+            h_max = parts[k].h;
             
     /* Get the cell width. */
-    if ( h_max < r_max )
-        h_max = r_max;
+    if ( h_cells < h_max )
+        h_cells = h_max;
     for ( k = 0 ; k < 3 ; k++ ) {
-        cdim[k] = ceil( dim[k] / h_max );
+        cdim[k] = ceil( dim[k] / h_cells );
         h[k] = dim[k] / cdim[k];
         ih[k] = 1.0 / h[k];
         }
@@ -1177,7 +1222,7 @@ void space_init ( struct space *s , double dim[3] , struct part *parts , int N ,
         }
         
     /* Store eveything in the space. */
-    s->r_min = r_min; s->r_max = r_max;
+    s->h_min = h_min; s->h_max = h_max;
     s->dim[0] = dim[0]; s->dim[1] = dim[1]; s->dim[2] = dim[2];
     s->periodic = periodic;
     s->parts = parts_new;
diff --git a/src/space.h b/src/space.h
index 61bf2e5503f52f0a21fc23eb5b533e897a224c0e..426d6ec04d0b592390f21f67190b09e2deaa26d5 100644
--- a/src/space.h
+++ b/src/space.h
@@ -52,7 +52,7 @@ struct space {
     double h[3], ih[3];
     
     /* The minimum and maximum cutoff radii. */
-    double r_min, r_max;
+    double h_min, h_max;
     
     /* Number of cells. */
     int nr_cells, tot_cells;
@@ -90,9 +90,10 @@ struct space {
 /* function prototypes. */
 struct cell *space_getcell ( struct space *s );
 struct task *space_gettask ( struct space *s );
+struct task *space_addtask ( struct space *s , int type , int subtype , int flags , int wait , struct cell *ci , struct cell *cj , struct task *unlock_tasks[] , int nr_unlock_tasks , struct cell *unlock_cells[] , int nr_unlock_cells );
 void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_max );
 void space_maketasks ( struct space *s , int do_sort );
-void space_map_cells ( struct space *s , void (*fun)( struct cell *c , void *data ) , void *data );
+void space_map_cells ( struct space *s , int full , void (*fun)( struct cell *c , void *data ) , void *data );
 void space_map_parts ( struct space *s , void (*fun)( struct part *p , struct cell *c , void *data ) , void *data );
 void space_recycle ( struct space *s , struct cell *c );
 
diff --git a/src/task.c b/src/task.c
index 4e1a7db6f716c2ce0ee9ebd10264adcff56b865c..6ba4f9b71db48ca495b79dbcc43ae3287a82ccfb 100644
--- a/src/task.c
+++ b/src/task.c
@@ -38,8 +38,7 @@
 
 
 /* Task type names. */
-const char *taskID_names[tid_count] = { "none" , "sort" , "self" , "pair" , "sub" };
-
+const char *taskID_names[task_type_count] = { "none" , "sort" , "self" , "pair" , "sub" , "ghost" };
 
 /* Error macro. */
 #define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
diff --git a/src/task.h b/src/task.h
index 7d5d3422b0c91c3eaee2becc04078cbc76dd1d5c..6ae474e116cc42eed4d99cb9485ca176ab3884ea 100644
--- a/src/task.h
+++ b/src/task.h
@@ -20,26 +20,36 @@
 
 /* Some constants. */
 #define task_maxwait                    3
-#define task_maxunlock                  39
+#define task_maxunlock                  40
 
 
-/* The different task IDs. */
-enum taskIDs {
-    tid_none = 0,
-    tid_sort,
-    tid_self,
-    tid_pair,
-    tid_sub,
-    tid_count
+/* The different task types. */
+enum task_types {
+    task_type_none = 0,
+    task_type_sort,
+    task_type_self,
+    task_type_pair,
+    task_type_sub,
+    task_type_ghost,
+    task_type_count
     };
     
 extern const char *taskID_names[];
     
+/* The different task sub-types. */
+enum task_subtypes {
+    task_subtype_none = 0,
+    task_subtype_density,
+    task_subtype_force,
+    task_subtype_count
+    };
+    
+extern const char *taskID_names[];
     
 /* Data of a task. */
 struct task {
 
-    int type, flags, wait, rank, done;
+    int type, subtype, flags, wait, rank, done;
     
     int nr_unlock_tasks;
     struct task *unlock_tasks[ task_maxunlock ];