diff --git a/CMakeLists.txt b/CMakeLists.txt
index c8d95552f1808ff52105afdd8e077fbbf3184ecd..bd6c307822e6ad078f265eea10ff4196da470e44 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -41,6 +41,13 @@ set(GPU_BIN "${TESTCASE}_gpu")
 
 set(RUNTIME_COMMON_FILES
     runtime/pairs.cpp
+    runtime/boundary_weights.cpp
+    runtime/copper_fcc_lattice.cpp
+    runtime/dem_sc_grid.cpp
+    runtime/read_from_file.cpp
+    runtime/stats.cpp
+    runtime/thermo.cpp
+    runtime/timing.cpp
     runtime/domain/block_forest.cpp
     runtime/domain/regular_6d_stencil.cpp)
 
diff --git a/runtime/boundary_weights.cpp b/runtime/boundary_weights.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..527789e5316d2f26f5285d8c8da1fc9beb4fa987
--- /dev/null
+++ b/runtime/boundary_weights.cpp
@@ -0,0 +1,74 @@
+#include <iostream>
+#include <string.h>
+#include <fstream>
+#include <sstream>
+//---
+#include "boundary_weights.hpp"
+#include "pairs.hpp"
+#include "pairs_common.hpp"
+
+// Always include last generated interfaces
+#include "interfaces/last_generated.hpp"
+
+#ifdef PAIRS_TARGET_CUDA
+int cuda_compute_boundary_weights(
+    real_t *position, int start, int end, int particle_capacity,
+    real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax);
+#endif
+
+namespace pairs {
+
+void compute_boundary_weights(
+    PairsRuntime *ps,
+    real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax,
+    long unsigned int *comp_weight, long unsigned int *comm_weight) {
+
+    const int particle_capacity = ps->getTrackedVariableAsInteger("particle_capacity");
+    const int nlocal = ps->getTrackedVariableAsInteger("nlocal");
+    const int nghost = ps->getTrackedVariableAsInteger("nghost");
+    auto position_prop = ps->getPropertyByName("position");
+
+    #ifndef PAIRS_TARGET_CUDA
+    real_t *position_ptr = static_cast<real_t *>(position_prop.getHostPointer());
+
+    *comp_weight = 0;
+    *comm_weight = 0;
+
+    for(int i = 0; i < nlocal; i++) {
+        real_t pos_x = pairs_host_interface::get_position(position_ptr, i, 0, particle_capacity);
+        real_t pos_y = pairs_host_interface::get_position(position_ptr, i, 1, particle_capacity);
+        real_t pos_z = pairs_host_interface::get_position(position_ptr, i, 2, particle_capacity);
+
+        if( pos_x > xmin && pos_x <= xmax &&
+            pos_y > ymin && pos_y <= ymax &&
+            pos_z > zmin && pos_z <= zmax) {
+                (*comp_weight)++;
+        }
+    }
+
+    for(int i = nlocal; i < nlocal + nghost; i++) {
+        real_t pos_x = pairs_host_interface::get_position(position_ptr, i, 0, particle_capacity);
+        real_t pos_y = pairs_host_interface::get_position(position_ptr, i, 1, particle_capacity);
+        real_t pos_z = pairs_host_interface::get_position(position_ptr, i, 2, particle_capacity);
+
+        if( pos_x > xmin && pos_x <= xmax &&
+            pos_y > ymin && pos_y <= ymax &&
+            pos_z > zmin && pos_z <= zmax) {
+                (*comm_weight)++;
+        }
+    }
+    std::cout << "comp_weight = " << (*comp_weight) << ", comm_weight = " << (*comm_weight) << std::endl;
+    #else
+    real_t *position_ptr = static_cast<real_t *>(position_prop.getDevicePointer());
+
+    ps->copyPropertyToDevice(position_prop, ReadOnly);
+
+    *comp_weight = cuda_compute_boundary_weights(
+        position_ptr, 0, nlocal, particle_capacity, xmin, xmax, ymin, ymax, zmin, zmax);
+
+    *comm_weight = cuda_compute_boundary_weights(
+        position_ptr, nlocal, nlocal + nghost, particle_capacity, xmin, xmax, ymin, ymax, zmin, zmax);
+    #endif
+}
+
+}
diff --git a/runtime/boundary_weights.hpp b/runtime/boundary_weights.hpp
index a679b59457802469b7d0657c1e5429ea03d8f8a9..c5d0e699c0d0e5ae7b32e0a97b35545717ed6c5f 100644
--- a/runtime/boundary_weights.hpp
+++ b/runtime/boundary_weights.hpp
@@ -1,10 +1,4 @@
-#include <iostream>
-#include <string.h>
-#include <fstream>
-#include <sstream>
-//---
 #include "pairs.hpp"
-#include "pairs_common.hpp"
 
 /*
 #define INTERFACE_DIR "interfaces/"
@@ -14,70 +8,13 @@
 #include INCLUDE_FILE(INTERFACE_FILE(INTERFACE_DIR, APPLICATION_REFERENCE, INTERFACE_EXT))
 */
 
-// Always include last generated interfaces
-#include "interfaces/last_generated.hpp"
-
 #pragma once
 
-#ifdef PAIRS_TARGET_CUDA
-int cuda_compute_boundary_weights(
-    real_t *position, int start, int end, int particle_capacity,
-    real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax);
-#endif
-
 namespace pairs {
 
 void compute_boundary_weights(
     PairsRuntime *ps,
     real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax,
-    walberla::uint_t *comp_weight, walberla::uint_t *comm_weight) {
-
-    const int particle_capacity = ps->getTrackedVariableAsInteger("particle_capacity");
-    const int nlocal = ps->getTrackedVariableAsInteger("nlocal");
-    const int nghost = ps->getTrackedVariableAsInteger("nghost");
-    auto position_prop = ps->getPropertyByName("position");
-
-    #ifndef PAIRS_TARGET_CUDA
-    real_t *position_ptr = static_cast<real_t *>(position_prop.getHostPointer());
-
-    *comp_weight = 0;
-    *comm_weight = 0;
-
-    for(int i = 0; i < nlocal; i++) {
-        real_t pos_x = pairs_host_interface::get_position(position_ptr, i, 0, particle_capacity);
-        real_t pos_y = pairs_host_interface::get_position(position_ptr, i, 1, particle_capacity);
-        real_t pos_z = pairs_host_interface::get_position(position_ptr, i, 2, particle_capacity);
-
-        if( pos_x > xmin && pos_x <= xmax &&
-            pos_y > ymin && pos_y <= ymax &&
-            pos_z > zmin && pos_z <= zmax) {
-                (*comp_weight)++;
-        }
-    }
-
-    for(int i = nlocal; i < nlocal + nghost; i++) {
-        real_t pos_x = pairs_host_interface::get_position(position_ptr, i, 0, particle_capacity);
-        real_t pos_y = pairs_host_interface::get_position(position_ptr, i, 1, particle_capacity);
-        real_t pos_z = pairs_host_interface::get_position(position_ptr, i, 2, particle_capacity);
-
-        if( pos_x > xmin && pos_x <= xmax &&
-            pos_y > ymin && pos_y <= ymax &&
-            pos_z > zmin && pos_z <= zmax) {
-                (*comm_weight)++;
-        }
-    }
-    std::cout << "comp_weight = " << (*comp_weight) << ", comm_weight = " << (*comm_weight) << std::endl;
-    #else
-    real_t *position_ptr = static_cast<real_t *>(position_prop.getDevicePointer());
-
-    ps->copyPropertyToDevice(position_prop, ReadOnly);
-
-    *comp_weight = cuda_compute_boundary_weights(
-        position_ptr, 0, nlocal, particle_capacity, xmin, xmax, ymin, ymax, zmin, zmax);
-
-    *comm_weight = cuda_compute_boundary_weights(
-        position_ptr, nlocal, nlocal + nghost, particle_capacity, xmin, xmax, ymin, ymax, zmin, zmax);
-    #endif
-}
+    long unsigned int *comp_weight, long unsigned int *comm_weight);
 
 }
diff --git a/runtime/copper_fcc_lattice.cpp b/runtime/copper_fcc_lattice.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..76c729c9e88830ffcf0274ea6fa434a8033b6d6b
--- /dev/null
+++ b/runtime/copper_fcc_lattice.cpp
@@ -0,0 +1,140 @@
+#include <iostream>
+#include <math.h>
+//---
+#include "copper_fcc_lattice.hpp"
+
+namespace pairs {
+
+double myrandom(int* seed) {
+    int k = (*seed) / IQ;
+    double ans;
+
+    *seed = IA * (*seed - k * IQ) - IR * k;
+    if(*seed < 0) *seed += IM;
+    ans = AM * (*seed);
+    return ans;
+}
+
+void random_reset(int *seed, int ibase, double *coord) {
+    int i;
+    char *str = (char *) &ibase;
+    int n = sizeof(int);
+    unsigned int hash = 0;
+
+    for (i = 0; i < n; i++) {
+        hash += str[i];
+        hash += (hash << 10);
+        hash ^= (hash >> 6);
+    }
+
+    str = (char *) coord;
+    n = 3 * sizeof(double);
+    for (i = 0; i < n; i++) {
+        hash += str[i];
+        hash += (hash << 10);
+        hash ^= (hash >> 6);
+    }
+
+    hash += (hash << 3);
+    hash ^= (hash >> 11);
+    hash += (hash << 15);
+
+    // keep 31 bits of unsigned int as new seed
+    // do not allow seed = 0, since will cause hang in gaussian()
+
+    *seed = hash & 0x7ffffff;
+    if (!(*seed)) *seed = 1;
+
+    // warm up the RNG
+
+    for (i = 0; i < 5; i++) myrandom(seed);
+    //save = 0;
+}
+
+double copper_fcc_lattice(
+    PairsRuntime *ps, int nx, int ny, int nz, double xprd, double yprd, double zprd,
+    double rho, int ntypes) {
+
+    auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+    auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
+    auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
+    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
+    auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
+    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+    double xlo = 0.0, xhi = xprd;
+    double ylo = 0.0, yhi = yprd;
+    double zlo = 0.0, zhi = zprd;
+    int natoms = 0;
+    //int natoms_expected = 4 * nx * ny * nz;
+
+    double alat = pow((4.0 / rho), (1.0 / 3.0));
+    int ilo = (int) (xlo / (0.5 * alat) - 1);
+    int ihi = (int) (xhi / (0.5 * alat) + 1);
+    int jlo = (int) (ylo / (0.5 * alat) - 1);
+    int jhi = (int) (yhi / (0.5 * alat) + 1);
+    int klo = (int) (zlo / (0.5 * alat) - 1);
+    int khi = (int) (zhi / (0.5 * alat) + 1);
+
+    ilo = MAX(ilo, 0);
+    ihi = MIN(ihi, 2 * nx - 1);
+    jlo = MAX(jlo, 0);
+    jhi = MIN(jhi, 2 * ny - 1);
+    klo = MAX(klo, 0);
+    khi = MIN(khi, 2 * nz - 1);
+
+    double xtmp, ytmp, ztmp, vxtmp, vytmp, vztmp;
+    int i, j, k, m, n;
+    int sx = 0; int sy = 0; int sz = 0;
+    int ox = 0; int oy = 0; int oz = 0;
+    int subboxdim = 8;
+
+    while(oz * subboxdim <= khi) {
+        k = oz * subboxdim + sz;
+        j = oy * subboxdim + sy;
+        i = ox * subboxdim + sx;
+
+        if(((i + j + k) % 2 == 0) &&
+            (i >= ilo) && (i <= ihi) &&
+            (j >= jlo) && (j <= jhi) &&
+            (k >= klo) && (k <= khi)) {
+
+            xtmp = 0.5 * alat * i;
+            ytmp = 0.5 * alat * j;
+            ztmp = 0.5 * alat * k;
+
+            if(ps->getDomainPartitioner()->isWithinSubdomain(xtmp, ytmp, ztmp)) {
+                n = k * (2 * ny) * (2 * nx) + j * (2 * nx) + i + 1;
+                for(m = 0; m < 5; m++) { myrandom(&n); }
+                vxtmp = myrandom(&n);
+                for(m = 0; m < 5; m++){ myrandom(&n); }
+                vytmp = myrandom(&n);
+                for(m = 0; m < 5; m++) { myrandom(&n); }
+                vztmp = myrandom(&n);
+
+                masses(natoms) = 1.0;
+                positions(natoms, 0) = xtmp;
+                positions(natoms, 1) = ytmp;
+                positions(natoms, 2) = ztmp;
+                velocities(natoms, 0) = vxtmp;
+                velocities(natoms, 1) = vytmp;
+                velocities(natoms, 2) = vztmp;
+                types(natoms) = rand() % ntypes;
+                flags(natoms) = 0;
+                shape(natoms) = 2; // point mass
+                natoms++;
+            }
+        }
+
+        sx++;
+
+        if(sx == subboxdim) { sx = 0; sy++; }
+        if(sy == subboxdim) { sy = 0; sz++; }
+        if(sz == subboxdim) { sz = 0; ox++; }
+        if(ox * subboxdim > ihi) { ox = 0; oy++; }
+        if(oy * subboxdim > jhi) { oy = 0; oz++; }
+    }
+
+    return natoms;
+}
+
+}
diff --git a/runtime/copper_fcc_lattice.hpp b/runtime/copper_fcc_lattice.hpp
index 254dd5d4a988e1d4f9e6f94c7369796895cb2040..39984e30deffc639ec5d327d25cf2357daab9357 100644
--- a/runtime/copper_fcc_lattice.hpp
+++ b/runtime/copper_fcc_lattice.hpp
@@ -1,12 +1,7 @@
-#include <iostream>
-#include <math.h>
-//---
 #include "pairs.hpp"
 
 #pragma once
 
-namespace pairs {
-
 /* Park/Miller RNG w/out MASKING, so as to be like f90s version */
 #define IA 16807
 #define IM 2147483647
@@ -15,133 +10,12 @@ namespace pairs {
 #define IR 2836
 #define MASK 123459876
 
-double myrandom(int* seed) {
-    int k= (*seed) / IQ;
-    double ans;
-
-    *seed = IA * (*seed - k * IQ) - IR * k;
-    if(*seed < 0) *seed += IM;
-    ans = AM * (*seed);
-    return ans;
-}
-
-void random_reset(int *seed, int ibase, double *coord) {
-    int i;
-    char *str = (char *) &ibase;
-    int n = sizeof(int);
-    unsigned int hash = 0;
-
-    for (i = 0; i < n; i++) {
-        hash += str[i];
-        hash += (hash << 10);
-        hash ^= (hash >> 6);
-    }
-
-    str = (char *) coord;
-    n = 3 * sizeof(double);
-    for (i = 0; i < n; i++) {
-        hash += str[i];
-        hash += (hash << 10);
-        hash ^= (hash >> 6);
-    }
-
-    hash += (hash << 3);
-    hash ^= (hash >> 11);
-    hash += (hash << 15);
-
-    // keep 31 bits of unsigned int as new seed
-    // do not allow seed = 0, since will cause hang in gaussian()
-
-    *seed = hash & 0x7ffffff;
-    if (!(*seed)) *seed = 1;
-
-    // warm up the RNG
-
-    for (i = 0; i < 5; i++) myrandom(seed);
-    //save = 0;
-}
-
-double copper_fcc_lattice(PairsRuntime *ps, int nx, int ny, int nz, double xprd, double yprd, double zprd, double rho, int ntypes) {
-    auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
-    auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
-    auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
-    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
-    auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
-    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
-    double xlo = 0.0, xhi = xprd;
-    double ylo = 0.0, yhi = yprd;
-    double zlo = 0.0, zhi = zprd;
-    int natoms = 0;
-    //int natoms_expected = 4 * nx * ny * nz;
-
-    double alat = pow((4.0 / rho), (1.0 / 3.0));
-    int ilo = (int) (xlo / (0.5 * alat) - 1);
-    int ihi = (int) (xhi / (0.5 * alat) + 1);
-    int jlo = (int) (ylo / (0.5 * alat) - 1);
-    int jhi = (int) (yhi / (0.5 * alat) + 1);
-    int klo = (int) (zlo / (0.5 * alat) - 1);
-    int khi = (int) (zhi / (0.5 * alat) + 1);
-
-    ilo = MAX(ilo, 0);
-    ihi = MIN(ihi, 2 * nx - 1);
-    jlo = MAX(jlo, 0);
-    jhi = MIN(jhi, 2 * ny - 1);
-    klo = MAX(klo, 0);
-    khi = MIN(khi, 2 * nz - 1);
-
-    double xtmp, ytmp, ztmp, vxtmp, vytmp, vztmp;
-    int i, j, k, m, n;
-    int sx = 0; int sy = 0; int sz = 0;
-    int ox = 0; int oy = 0; int oz = 0;
-    int subboxdim = 8;
-
-    while(oz * subboxdim <= khi) {
-        k = oz * subboxdim + sz;
-        j = oy * subboxdim + sy;
-        i = ox * subboxdim + sx;
-
-        if(((i + j + k) % 2 == 0) &&
-            (i >= ilo) && (i <= ihi) &&
-            (j >= jlo) && (j <= jhi) &&
-            (k >= klo) && (k <= khi)) {
-
-            xtmp = 0.5 * alat * i;
-            ytmp = 0.5 * alat * j;
-            ztmp = 0.5 * alat * k;
-
-            if(ps->getDomainPartitioner()->isWithinSubdomain(xtmp, ytmp, ztmp)) {
-                n = k * (2 * ny) * (2 * nx) + j * (2 * nx) + i + 1;
-                for(m = 0; m < 5; m++) { myrandom(&n); }
-                vxtmp = myrandom(&n);
-                for(m = 0; m < 5; m++){ myrandom(&n); }
-                vytmp = myrandom(&n);
-                for(m = 0; m < 5; m++) { myrandom(&n); }
-                vztmp = myrandom(&n);
-
-                masses(natoms) = 1.0;
-                positions(natoms, 0) = xtmp;
-                positions(natoms, 1) = ytmp;
-                positions(natoms, 2) = ztmp;
-                velocities(natoms, 0) = vxtmp;
-                velocities(natoms, 1) = vytmp;
-                velocities(natoms, 2) = vztmp;
-                types(natoms) = rand() % ntypes;
-                flags(natoms) = 0;
-                shape(natoms) = 2; // point mass
-                natoms++;
-            }
-        }
-
-        sx++;
-
-        if(sx == subboxdim) { sx = 0; sy++; }
-        if(sy == subboxdim) { sy = 0; sz++; }
-        if(sz == subboxdim) { sz = 0; ox++; }
-        if(ox * subboxdim > ihi) { ox = 0; oy++; }
-        if(oy * subboxdim > jhi) { oy = 0; oz++; }
-    }
+namespace pairs {
 
-    return natoms;
-}
+double myrandom(int* seed);
+void random_reset(int *seed, int ibase, double *coord);
+double copper_fcc_lattice(
+    PairsRuntime *ps, int nx, int ny, int nz, double xprd, double yprd, double zprd,
+    double rho, int ntypes);
 
 }
diff --git a/runtime/dem_sc_grid.cpp b/runtime/dem_sc_grid.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..66aa5ac59a0803b0a8a2294160ec034dce48cc5a
--- /dev/null
+++ b/runtime/dem_sc_grid.cpp
@@ -0,0 +1,140 @@
+#include <iostream>
+#include <math.h>
+#include <random>
+//---
+#include "dem_sc_grid.hpp"
+#include "pairs.hpp"
+#include "pairs_common.hpp"
+
+namespace pairs {
+
+namespace internal {
+
+static std::mt19937 generator; // static std::mt19937_64 generator;
+
+std::mt19937 & get_generator() {
+    // std::mt19937_64
+    return generator;
+}
+
+}
+
+bool point_within_aabb(double point[], double aabb[]) {
+    return point[0] >= aabb[0] && point[0] < aabb[3] &&
+           point[1] >= aabb[1] && point[1] < aabb[4] &&
+           point[2] >= aabb[2] && point[2] < aabb[5];
+}
+
+int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double spacing, double diameter, double min_diameter, double max_diameter, double initial_velocity, double particle_density, int ntypes) {
+    auto uid = ps->getAsIntegerProperty(ps->getPropertyByName("uid"));
+    auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+    auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
+    auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
+    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
+    auto radius = ps->getAsFloatProperty(ps->getPropertyByName("radius"));
+    auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
+    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+    int last_uid = 1;
+    int nparticles = 0;
+
+    const double xmin = 0.0;
+    const double ymin = 0.0;
+    const double zmin = diameter;
+
+    double gen_domain[] = {xmin, ymin, zmin, xmax, ymax, zmax};
+    double ref_point[] = {spacing * 0.5, spacing * 0.5, spacing * 0.5};
+    double sc_xmin = xmin - ref_point[0];
+    double sc_ymin = ymin - ref_point[1];
+    double sc_zmin = zmin - ref_point[2];
+
+    int iret = (int)(ceil(sc_xmin / spacing));
+    int jret = (int)(ceil(sc_ymin / spacing));
+    int kret = (int)(ceil(sc_zmin / spacing));
+
+    int i = iret;
+    int j = jret;
+    int k = kret;
+
+    double point[3];
+    point[0] = ref_point[0] + i * spacing;
+    point[1] = ref_point[1] + j * spacing;
+    point[2] = ref_point[2] + k * spacing;
+
+    while(point_within_aabb(point, gen_domain)) {
+        int particle_uid = last_uid;
+        auto pdiam = realRandom<real_t>(min_diameter, max_diameter);
+
+        if(ps->getDomainPartitioner()->isWithinSubdomain(point[0], point[1], point[2])) {
+            real_t rad = pdiam * 0.5;
+            uid(nparticles) = particle_uid;
+            radius(nparticles) = rad;
+            masses(nparticles) = ((4.0 / 3.0) * M_PI) * rad * rad * rad * particle_density;
+            positions(nparticles, 0) = point[0];
+            positions(nparticles, 1) = point[1];
+            positions(nparticles, 2) = point[2];
+            velocities(nparticles, 0) = 0.1 * realRandom<real_t>(-initial_velocity, initial_velocity);
+            velocities(nparticles, 1) = 0.1 * realRandom<real_t>(-initial_velocity, initial_velocity);
+            velocities(nparticles, 2) = -initial_velocity;
+            types(nparticles) = rand() % ntypes;
+            flags(nparticles) = 0;
+            shape(nparticles) = 0; // sphere
+
+            /*
+            std::cout << uid(nparticles) << "," << types(nparticles) << "," << masses(nparticles) << "," << radius(nparticles) << ","
+                      << positions(nparticles, 0) << "," << positions(nparticles, 1) << "," << positions(nparticles, 2) << ","
+                      << velocities(nparticles, 0) << "," << velocities(nparticles, 1) << "," << velocities(nparticles, 2) << ","
+                      << flags(nparticles) << std::endl;
+            */
+
+            nparticles++;
+        }
+
+        ++i;
+        point[0] = ref_point[0] + i * spacing;
+        point[1] = ref_point[1] + j * spacing;
+        point[2] = ref_point[2] + k * spacing;
+
+        if(!point_within_aabb(point, gen_domain)) {
+            i = iret;
+            j++;
+            point[0] = ref_point[0] + i * spacing;
+            point[1] = ref_point[1] + j * spacing;
+            point[2] = ref_point[2] + k * spacing;
+
+            if(!point_within_aabb(point, gen_domain)) {
+                j = jret;
+                k++;
+                point[0] = ref_point[0] + i * spacing;
+                point[1] = ref_point[1] + j * spacing;
+                point[2] = ref_point[2] + k * spacing;
+
+                if(!point_within_aabb(point, gen_domain)) {
+                    break;
+                }
+            }
+        }
+
+        last_uid++;
+    }
+
+    int global_nparticles = nparticles;
+    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+        MPI_Allreduce(&nparticles, &global_nparticles, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+    }
+
+    if(ps->getDomainPartitioner()->getRank() == 0) {
+        std::cout << "DEM Simple-Cubic Grid" << std::endl;
+        std::cout << "Domain size: <" << xmax << ", " << ymax << ", " << zmax << ">" << std::endl;
+        std::cout << "Spacing: " << spacing << std::endl;
+        std::cout << "Diameter: " << diameter
+                  << " (min = " << min_diameter << ", max = " << max_diameter << ")" << std::endl;
+        std::cout << "Initial velocity: " << initial_velocity << std::endl;
+        std::cout << "Particle density: " << particle_density << std::endl;
+        std::cout << "Number of types: " << ntypes << std::endl;
+        std::cout << "Number of particles: " << global_nparticles << std::endl;
+    }
+
+    return nparticles;
+}
+
+}
diff --git a/runtime/dem_sc_grid.hpp b/runtime/dem_sc_grid.hpp
index 8d4848e7c9914501ed2da734b99a696f87edbf02..a293433095d285441abacbdbad09533148133c0f 100644
--- a/runtime/dem_sc_grid.hpp
+++ b/runtime/dem_sc_grid.hpp
@@ -1,8 +1,8 @@
-#include <iostream>
 #include <math.h>
 #include <random>
 //---
 #include "pairs.hpp"
+#include "pairs_common.hpp"
 
 #pragma once
 
@@ -10,16 +10,11 @@ namespace pairs {
 
 namespace internal {
 
-static std::mt19937 generator; // static std::mt19937_64 generator;
+std::mt19937 & get_generator();
 
-std::mt19937 & get_generator() {
-    // std::mt19937_64
-    return generator;
 }
 
-}
-
-template< typename REAL_TYPE = real_t>
+template<typename REAL_TYPE = real_t>
 REAL_TYPE realRandom(
     const REAL_TYPE min = REAL_TYPE(0),
     const REAL_TYPE max = REAL_TYPE(1),
@@ -41,134 +36,23 @@ REAL_TYPE realRandom(
    return value;
 }
 
+template<typename REAL_TYPE>
+class RealRandom {
+public:
+    RealRandom(const std::mt19937::result_type& seed = std::mt19937::result_type()) {
+        generator_.seed(seed);
+    }
 
+    REAL_TYPE operator()(const REAL_TYPE min = REAL_TYPE(0), const REAL_TYPE max = REAL_TYPE(1)) {
+        return realRandom(min, max, generator_);
+    }
 
-template<typename REAL_TYPE> class RealRandom {
-public:
-   RealRandom(const std::mt19937::result_type& seed = std::mt19937::result_type()) { generator_.seed(seed); }
-   REAL_TYPE operator()(const REAL_TYPE min = REAL_TYPE(0), const REAL_TYPE max = REAL_TYPE(1) ) {
-      return realRandom(min, max, generator_);
-   }
 private:
    std::mt19937 generator_;
 };
 
-bool point_within_aabb(double point[], double aabb[]) {
-    return point[0] >= aabb[0] && point[0] < aabb[3] &&
-           point[1] >= aabb[1] && point[1] < aabb[4] &&
-           point[2] >= aabb[2] && point[2] < aabb[5];
-}
-
-int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double spacing, double diameter, double min_diameter, double max_diameter, double initial_velocity, double particle_density, int ntypes) {
-    auto uid = ps->getAsIntegerProperty(ps->getPropertyByName("uid"));
-    auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
-    auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
-    auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
-    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
-    auto radius = ps->getAsFloatProperty(ps->getPropertyByName("radius"));
-    auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
-    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
-    int last_uid = 1;
-    int nparticles = 0;
-
-    const double xmin = 0.0;
-    const double ymin = 0.0;
-    const double zmin = diameter;
-
-    double gen_domain[] = {xmin, ymin, zmin, xmax, ymax, zmax};
-    double ref_point[] = {spacing * 0.5, spacing * 0.5, spacing * 0.5};
-    double sc_xmin = xmin - ref_point[0];
-    double sc_ymin = ymin - ref_point[1];
-    double sc_zmin = zmin - ref_point[2];
-
-    int iret = (int)(ceil(sc_xmin / spacing));
-    int jret = (int)(ceil(sc_ymin / spacing));
-    int kret = (int)(ceil(sc_zmin / spacing));
-
-    int i = iret;
-    int j = jret;
-    int k = kret;
-
-    double point[3];
-    point[0] = ref_point[0] + i * spacing;
-    point[1] = ref_point[1] + j * spacing;
-    point[2] = ref_point[2] + k * spacing;
-
-    while(point_within_aabb(point, gen_domain)) {
-        int particle_uid = last_uid;
-        auto diameter = realRandom<real_t>(min_diameter, max_diameter);
-
-        if(ps->getDomainPartitioner()->isWithinSubdomain(point[0], point[1], point[2])) {
-            real_t rad = diameter * 0.5;
-            uid(nparticles) = particle_uid;
-            radius(nparticles) = rad;
-            masses(nparticles) = ((4.0 / 3.0) * M_PI) * rad * rad * rad * particle_density;
-            positions(nparticles, 0) = point[0];
-            positions(nparticles, 1) = point[1];
-            positions(nparticles, 2) = point[2];
-            velocities(nparticles, 0) = 0.1 * realRandom<real_t>(-initial_velocity, initial_velocity);
-            velocities(nparticles, 1) = 0.1 * realRandom<real_t>(-initial_velocity, initial_velocity);
-            velocities(nparticles, 2) = -initial_velocity;
-            types(nparticles) = rand() % ntypes;
-            flags(nparticles) = 0;
-            shape(nparticles) = 0; // sphere
-
-            /*
-            std::cout << uid(nparticles) << "," << types(nparticles) << "," << masses(nparticles) << "," << radius(nparticles) << ","
-                      << positions(nparticles, 0) << "," << positions(nparticles, 1) << "," << positions(nparticles, 2) << ","
-                      << velocities(nparticles, 0) << "," << velocities(nparticles, 1) << "," << velocities(nparticles, 2) << ","
-                      << flags(nparticles) << std::endl;
-            */
-
-            nparticles++;
-        }
-
-        ++i;
-        point[0] = ref_point[0] + i * spacing;
-        point[1] = ref_point[1] + j * spacing;
-        point[2] = ref_point[2] + k * spacing;
-
-        if(!point_within_aabb(point, gen_domain)) {
-            i = iret;
-            j++;
-            point[0] = ref_point[0] + i * spacing;
-            point[1] = ref_point[1] + j * spacing;
-            point[2] = ref_point[2] + k * spacing;
-
-            if(!point_within_aabb(point, gen_domain)) {
-                j = jret;
-                k++;
-                point[0] = ref_point[0] + i * spacing;
-                point[1] = ref_point[1] + j * spacing;
-                point[2] = ref_point[2] + k * spacing;
-
-                if(!point_within_aabb(point, gen_domain)) {
-                    break;
-                }
-            }
-        }
-
-        last_uid++;
-    }
-
-    int global_nparticles = nparticles;
-    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
-        MPI_Allreduce(&nparticles, &global_nparticles, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-    }
+bool point_within_aabb(double point[], double aabb[]);
 
-    if(ps->getDomainPartitioner()->getRank() == 0) {
-        std::cout << "DEM Simple-Cubic Grid" << std::endl;
-        std::cout << "Domain size: <" << xmax << ", " << ymax << ", " << zmax << ">" << std::endl;
-        std::cout << "Spacing: " << spacing << std::endl;
-        std::cout << "Diameter: " << diameter
-                  << " (min = " << min_diameter << ", max = " << max_diameter << ")" << std::endl;
-        std::cout << "Initial velocity: " << initial_velocity << std::endl;
-        std::cout << "Particle density: " << particle_density << std::endl;
-        std::cout << "Number of types: " << ntypes << std::endl;
-        std::cout << "Number of particles: " << global_nparticles << std::endl;
-    }
-
-    return nparticles;
-}
+int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double spacing, double diameter, double min_diameter, double max_diameter, double initial_velocity, double particle_density, int ntypes);
 
 }
diff --git a/runtime/device_flags.hpp b/runtime/device_flags.hpp
index 4b5085fb6b61f8df1c3e1e129541266765745fdf..089e32f3122dff09e56e01120ca33816844f8030 100644
--- a/runtime/device_flags.hpp
+++ b/runtime/device_flags.hpp
@@ -14,7 +14,7 @@ private:
     static const int narrays_per_flag = 64;
 public:
     DeviceFlags(int narrays_) : narrays(narrays_) {
-        nflags = std::ceil((double) narrays_ / (double) narrays_per_flag);
+        nflags = static_cast<int>(std::ceil((double) narrays_ / (double) narrays_per_flag));
         hflags = new unsigned long long int[nflags];
         dflags = new unsigned long long int[nflags];
 
diff --git a/runtime/devices/device.hpp b/runtime/devices/device.hpp
index 107b70ee91512ed9ccd336be2168e0b75ed5eab8..e340f12ba4276ac1b7d4e5d67fda245acb115d18 100644
--- a/runtime/devices/device.hpp
+++ b/runtime/devices/device.hpp
@@ -103,7 +103,7 @@ __device__ int atomic_add_resize_check(int *addr, int val, int *resize, int capa
 }
 #else
 inline int atomic_add(int *addr, int val) { return host_atomic_add(addr, val); }
-inline int atomic_add(real_t *addr, real_t val) { return host_atomic_add(addr, val); }
+inline real_t atomic_add(real_t *addr, real_t val) { return host_atomic_add(addr, val); }
 inline int atomic_add_resize_check(int *addr, int val, int *resize, int capacity) {
     return host_atomic_add_resize_check(addr, val, resize, capacity);
 }
diff --git a/runtime/domain/block_forest.hpp b/runtime/domain/block_forest.hpp
index 25712c73a5651b76c6242b300990f3581f057b0a..15aac5b59519f87d8e0ae85a6633cf88ecdfd910 100644
--- a/runtime/domain/block_forest.hpp
+++ b/runtime/domain/block_forest.hpp
@@ -22,13 +22,13 @@ class PairsRuntime;
 
 class BlockForest : public DomainPartitioner {
 private:
-    PairsRuntime *ps;
     std::shared_ptr<walberla::BlockForest> forest;
     std::shared_ptr<walberla::blockforest::InfoCollection> info;
     std::vector<int> ranks;
     std::vector<int> naabbs;
     std::vector<int> aabb_offsets;
     std::vector<double> aabbs;
+    PairsRuntime *ps;
     real_t *subdom;
     int world_size, rank, nranks, total_aabbs;
     bool balance_workload;
@@ -37,8 +37,7 @@ public:
     BlockForest(
         PairsRuntime *ps_,
         real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) :
-        ps(ps_),
-        DomainPartitioner(xmin, xmax, ymin, ymax, zmin, zmax) {
+        DomainPartitioner(xmin, xmax, ymin, ymax, zmin, zmax), ps(ps_) {
 
         subdom = new real_t[ndims * 2];
     }
diff --git a/runtime/pairs.cpp b/runtime/pairs.cpp
index 7cad7d1fc0f58bbb0717d621596c9905a78808ba..958f781e983fa4e26d32356398935b608e99264b 100644
--- a/runtime/pairs.cpp
+++ b/runtime/pairs.cpp
@@ -32,13 +32,11 @@ void PairsRuntime::initDomain(
 }
 
 void PairsRuntime::addArray(Array array) {
-    int id = array.getId();
-    auto a = std::find_if(
-        arrays.begin(),
-        arrays.end(),
-        [id](Array _a) { return _a.getId() == id; });
+    PAIRS_ASSERT(
+        std::find_if(arrays.begin(), arrays.end(), [array](Array _a) {
+            return _a.getId() == array.getId();
+        }) == std::end(arrays));
 
-    PAIRS_ASSERT(a == std::end(arrays));
     arrays.push_back(array);
 }
 
@@ -73,13 +71,11 @@ Array &PairsRuntime::getArrayByHostPointer(const void *h_ptr) {
 }
 
 void PairsRuntime::addProperty(Property prop) {
-    int id = prop.getId();
-    auto p = std::find_if(
-        properties.begin(),
-        properties.end(),
-        [id](Property _p) { return _p.getId() == id; });
+    PAIRS_ASSERT(
+        std::find_if(properties.begin(), properties.end(), [prop](Property _p) {
+            return _p.getId() == prop.getId();
+        }) == std::end(properties));
 
-    PAIRS_ASSERT(p == std::end(properties));
     properties.push_back(prop);
 }
 
@@ -104,13 +100,14 @@ Property &PairsRuntime::getPropertyByName(std::string name) {
 }
 
 void PairsRuntime::addContactProperty(ContactProperty contact_prop) {
-    int id = contact_prop.getId();
-    auto cp = std::find_if(
-        contact_properties.begin(),
-        contact_properties.end(),
-        [id](ContactProperty _cp) { return _cp.getId() == id; });
+    PAIRS_ASSERT(
+        std::find_if(
+            contact_properties.begin(),
+            contact_properties.end(),
+            [contact_prop](ContactProperty _cp) {
+                return _cp.getId() == contact_prop.getId();
+            }) == std::end(contact_properties));
 
-    PAIRS_ASSERT(cp == std::end(contact_properties));
     contact_properties.push_back(contact_prop);
 }
 
@@ -135,13 +132,14 @@ ContactProperty &PairsRuntime::getContactPropertyByName(std::string name) {
 }
 
 void PairsRuntime::addFeatureProperty(FeatureProperty feature_prop) {
-    int id = feature_prop.getId();
-    auto fp = std::find_if(
-        feature_properties.begin(),
-        feature_properties.end(),
-        [id](FeatureProperty _fp) { return _fp.getId() == id; });
+    PAIRS_ASSERT(
+        std::find_if(
+            feature_properties.begin(),
+            feature_properties.end(),
+            [feature_prop](FeatureProperty _fp) {
+                return _fp.getId() == feature_prop.getId();
+            }) == std::end(feature_properties));
 
-    PAIRS_ASSERT(fp == std::end(feature_properties));
     feature_properties.push_back(feature_prop);
 }
 
@@ -170,7 +168,7 @@ void PairsRuntime::copyArraySliceToDevice(
         if(action == Ignore || !array_flags->isDeviceFlagSet(array_id)) {
             if(!array.isStatic()) {
                 PAIRS_DEBUG(
-                    "Copying array %s to device (offset=%d, n=%d)\n",
+                    "Copying array %s to device (offset=%lu, n=%lu)\n",
                     array.getName().c_str(), offset, size);
 
                 pairs::copy_slice_to_device(
@@ -192,10 +190,17 @@ void PairsRuntime::copyArrayToDevice(Array &array, action_t action, size_t size)
     if(action == Ignore || action == WriteAfterRead || action == ReadOnly) {
         if(action == Ignore || !array_flags->isDeviceFlagSet(array_id)) {
             if(array.isStatic()) {
-                PAIRS_DEBUG("Copying static array %s to device (n=%d)\n", array.getName().c_str(), size);
-                pairs::copy_static_symbol_to_device(array.getHostPointer(), array.getDevicePointer(), size);
+                PAIRS_DEBUG(
+                    "Copying static array %s to device (n=%lu)\n",
+                    array.getName().c_str(), size);
+
+                pairs::copy_static_symbol_to_device(
+                    array.getHostPointer(), array.getDevicePointer(), size);
             } else {
-                PAIRS_DEBUG("Copying array %s to device (n=%d)\n", array.getName().c_str(), size);
+                PAIRS_DEBUG(
+                    "Copying array %s to device (n=%lu)\n",
+                    array.getName().c_str(), size);
+
                 pairs::copy_to_device(array.getHostPointer(), array.getDevicePointer(), size);
             }
         }
@@ -215,7 +220,7 @@ void PairsRuntime::copyArraySliceToHost(Array &array, action_t action, size_t of
         if(action == Ignore || !array_flags->isHostFlagSet(array_id)) {
             if(!array.isStatic()) {
                 PAIRS_DEBUG(
-                    "Copying array %s to host (offset=%d, n=%d)\n",
+                    "Copying array %s to host (offset=%lu, n=%lu)\n",
                     array.getName().c_str(), offset, size);
 
                 pairs::copy_slice_to_host(
@@ -237,10 +242,13 @@ void PairsRuntime::copyArrayToHost(Array &array, action_t action, size_t size) {
     if(action == Ignore || action == WriteAfterRead || action == ReadOnly) {
         if(action == Ignore || !array_flags->isHostFlagSet(array_id)) {
             if(array.isStatic()) {
-                PAIRS_DEBUG("Copying static array %s to host (n=%d)\n", array.getName().c_str(), size);
-                pairs::copy_static_symbol_to_host(array.getDevicePointer(), array.getHostPointer(), size);
+                PAIRS_DEBUG(
+                    "Copying static array %s to host (n=%lu)\n", array.getName().c_str(), size);
+
+                pairs::copy_static_symbol_to_host(
+                    array.getDevicePointer(), array.getHostPointer(), size);
             } else {
-                PAIRS_DEBUG("Copying array %s to host (n=%d)\n", array.getName().c_str(), size);
+                PAIRS_DEBUG("Copying array %s to host (n=%lu)\n", array.getName().c_str(), size);
                 pairs::copy_to_host(array.getDevicePointer(), array.getHostPointer(), size);
             }
         }
@@ -258,7 +266,7 @@ void PairsRuntime::copyPropertyToDevice(Property &prop, action_t action, size_t
 
     if(action == Ignore || action == WriteAfterRead || action == ReadOnly) {
         if(action == Ignore || !prop_flags->isDeviceFlagSet(prop_id)) {
-            PAIRS_DEBUG("Copying property %s to device (n=%d)\n", prop.getName().c_str(), size);
+            PAIRS_DEBUG("Copying property %s to device (n=%lu)\n", prop.getName().c_str(), size);
             pairs::copy_to_device(prop.getHostPointer(), prop.getDevicePointer(), size);
         }
     }
@@ -275,7 +283,7 @@ void PairsRuntime::copyPropertyToHost(Property &prop, action_t action, size_t si
 
     if(action == Ignore || action == WriteAfterRead || action == ReadOnly) {
         if(action == Ignore || !prop_flags->isHostFlagSet(prop_id)) {
-            PAIRS_DEBUG("Copying property %s to host (n=%d)\n", prop.getName().c_str(), size);
+            PAIRS_DEBUG("Copying property %s to host (n=%lu)\n", prop.getName().c_str(), size);
             pairs::copy_to_host(prop.getDevicePointer(), prop.getHostPointer(), size);
         }
     }
@@ -294,8 +302,12 @@ void PairsRuntime::copyContactPropertyToDevice(
 
     if(action == Ignore || action == WriteAfterRead || action == ReadOnly) {
         if(action == Ignore || !contact_prop_flags->isDeviceFlagSet(prop_id)) {
-            PAIRS_DEBUG("Copying contact property %s to device (n=%d)\n", contact_prop.getName().c_str(), size);
-            pairs::copy_to_device(contact_prop.getHostPointer(), contact_prop.getDevicePointer(), size);
+            PAIRS_DEBUG("Copying contact property %s to device (n=%lu)\n",
+                contact_prop.getName().c_str(), size);
+
+            pairs::copy_to_device(
+                contact_prop.getHostPointer(), contact_prop.getDevicePointer(), size);
+
             contact_prop_flags->setDeviceFlag(prop_id);
         }
     }
@@ -312,8 +324,12 @@ void PairsRuntime::copyContactPropertyToHost(
 
     if(action == Ignore || action == WriteAfterRead || action == ReadOnly) {
         if(!contact_prop_flags->isHostFlagSet(contact_prop.getId())) {
-            PAIRS_DEBUG("Copying contact property %s to host (n=%d)\n", contact_prop.getName().c_str(), size);
-            pairs::copy_to_host(contact_prop.getDevicePointer(), contact_prop.getHostPointer(), size);
+            PAIRS_DEBUG("Copying contact property %s to host (n=%lu)\n",
+                contact_prop.getName().c_str(), size);
+
+            pairs::copy_to_host(
+                contact_prop.getDevicePointer(), contact_prop.getHostPointer(), size);
+
             contact_prop_flags->setHostFlag(prop_id);
         }
     }
@@ -325,8 +341,12 @@ void PairsRuntime::copyContactPropertyToHost(
 
 void PairsRuntime::copyFeaturePropertyToDevice(FeatureProperty &feature_prop) {
     const size_t n = feature_prop.getArraySize();
-    PAIRS_DEBUG("Copying feature property %s to device (n=%d)\n", feature_prop.getName().c_str(), n);
-    pairs::copy_static_symbol_to_device(feature_prop.getHostPointer(), feature_prop.getDevicePointer(), n);
+
+    PAIRS_DEBUG("Copying feature property %s to device (n=%lu)\n",
+        feature_prop.getName().c_str(), n);
+
+    pairs::copy_static_symbol_to_device(
+        feature_prop.getHostPointer(), feature_prop.getDevicePointer(), n);
 }
 
 void PairsRuntime::communicateSizes(int dim, const int *send_sizes, int *recv_sizes) {
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index bb8375f79de2b95fbb20fc0bcdee1af090f41d18..d7716131f2f3d9942607d90472a0534fe1e9a538 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -68,12 +68,12 @@ public:
     }
 
     void trackVariable(std::string variable_name, void *ptr) {
-        auto v = std::find_if(
-            tracked_variables.begin(),
-            tracked_variables.end(),
-            [variable_name](TrackedVariable _v) { return _v.getName() == variable_name; });
+        PAIRS_ASSERT(
+            std::find_if(tracked_variables.begin(), tracked_variables.end(),
+            [variable_name](TrackedVariable _v) {
+                return _v.getName() == variable_name;
+            }) == std::end(tracked_variables));
 
-        PAIRS_ASSERT(v == std::end(tracked_variables));
         tracked_variables.push_back(TrackedVariable(variable_name, ptr)); 
     }
 
@@ -83,7 +83,7 @@ public:
             tracked_variables.end(),
             [variable_name](TrackedVariable _v) { return _v.getName() == variable_name; });
 
-        PAIRS_ASSERT(v == std::end(tracked_variables));
+        PAIRS_ASSERT(v != std::end(tracked_variables));
         return *v;
     }
 
@@ -92,7 +92,7 @@ public:
         *(static_cast<int *>(tv.getPointer())) = value;
     }
 
-    const int getTrackedVariableAsInteger(std::string variable_name) {
+    int getTrackedVariableAsInteger(std::string variable_name) {
         auto& tv = getTrackedVariable(variable_name);
         return *(static_cast<int *>(tv.getPointer()));
     }
@@ -371,7 +371,7 @@ void PairsRuntime::addStaticArray(array_t id, std::string name, T_ptr *h_ptr, T_
 template<typename T_ptr>
 void PairsRuntime::reallocArray(array_t id, T_ptr **h_ptr, std::nullptr_t, size_t size) {
     // This should be a pointer (and not a reference) in order to be modified
-    auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
+    auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array _a) { return _a.getId() == id; });
     PAIRS_ASSERT(a != std::end(arrays));
     PAIRS_ASSERT(size > 0);
 
@@ -386,7 +386,7 @@ void PairsRuntime::reallocArray(array_t id, T_ptr **h_ptr, std::nullptr_t, size_
 template<typename T_ptr>
 void PairsRuntime::reallocArray(array_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t size) {
     // This should be a pointer (and not a reference) in order to be modified
-    auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
+    auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array _a) { return _a.getId() == id; });
     PAIRS_ASSERT(a != std::end(arrays));
     PAIRS_ASSERT(size > 0);
 
diff --git a/runtime/property.hpp b/runtime/property.hpp
index 84a2bc38151ee8ada0a0ebcf2283f8d5d5303bc5..13ef0460b9295e4d300c5669aafddf46fac995a3 100644
--- a/runtime/property.hpp
+++ b/runtime/property.hpp
@@ -11,8 +11,8 @@ protected:
     void *h_ptr, *d_ptr;
     PropertyType type;
     layout_t layout;
-    size_t sx, sy;
     int vol;
+    size_t sx, sy;
 
 public:
     Property(
diff --git a/runtime/read_from_file.cpp b/runtime/read_from_file.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..73058071a6f320b3baabdcf30c3324005ad810c3
--- /dev/null
+++ b/runtime/read_from_file.cpp
@@ -0,0 +1,168 @@
+#include <iostream>
+#include <string.h>
+#include <fstream>
+#include <sstream>
+//---
+#include "pairs.hpp"
+#include "pairs_common.hpp"
+
+namespace pairs {
+
+void read_grid_data(PairsRuntime *ps, const char *filename, real_t *grid_buffer) {
+    std::ifstream in_file(filename, std::ifstream::in);
+    std::string line;
+
+    if(in_file.is_open()) {
+        std::getline(in_file, line);
+        std::stringstream line_stream(line);
+        std::string in0;
+        int i = 0;
+
+        while(std::getline(line_stream, in0, ',')) {
+            //PAIRS_ASSERT(i < ndims * 2);
+            grid_buffer[i] = std::stod(in0);
+            i++;
+        }
+
+        in_file.close();
+    }
+}
+
+size_t read_particle_data(
+    PairsRuntime *ps, const char *filename, const property_t properties[],
+    size_t nprops, int shape_id, int start) {
+
+    std::ifstream in_file(filename, std::ifstream::in);
+    std::string line;
+    auto shape_ptr = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+    int n = start;
+
+    if(in_file.is_open()) {
+        //std::getline(in_file, line);
+        while(std::getline(in_file, line)) {
+            std::stringstream line_stream(line);
+            std::string in0;
+            int within_domain = 1;
+            int i = 0;
+            int flags = 0;
+
+            while(std::getline(line_stream, in0, ',')) {
+                property_t p_id = properties[i];
+                auto prop = ps->getProperty(p_id);
+                auto prop_type = prop.getType();
+
+                if(prop_type == Prop_Vector) {
+                    auto vector_ptr = ps->getAsVectorProperty(prop);
+                    std::string in1, in2;
+                    std::getline(line_stream, in1, ',');
+                    std::getline(line_stream, in2, ',');
+                    real_t x = std::stod(in0);
+                    real_t y = std::stod(in1);
+                    real_t z = std::stod(in2);
+                    vector_ptr(n, 0) = x;
+                    vector_ptr(n, 1) = y;
+                    vector_ptr(n, 2) = z;
+
+                    if(prop.getName() == "position") {
+                        within_domain = ps->getDomainPartitioner()->isWithinSubdomain(x, y, z);
+                    }
+                } else if(prop_type == Prop_Matrix) {
+                    auto matrix_ptr = ps->getAsMatrixProperty(prop);
+                    constexpr int nelems = 9;
+                    std::string in_buf;
+
+                    matrix_ptr(n, 0) = std::stod(in0);
+                    for(int e = 1; e < nelems; e++) {
+                        std::getline(line_stream, in_buf, ',');
+                        matrix_ptr(n, e) = std::stod(in_buf);
+                    }
+                } else if(prop_type == Prop_Quaternion) {
+                    auto quat_ptr = ps->getAsQuaternionProperty(prop);
+                    constexpr int nelems = 4;
+                    std::string in_buf;
+
+                    quat_ptr(n, 0) = std::stod(in0);
+                    for(int e = 1; e < nelems; e++) {
+                        std::getline(line_stream, in_buf, ',');
+                        quat_ptr(n, e) = std::stod(in_buf);
+                    }
+                } else if(prop_type == Prop_Integer) {
+                    auto int_ptr = ps->getAsIntegerProperty(prop);
+                    int_ptr(n) = std::stoi(in0);
+
+                    if(prop.getName() == "flags") {
+                        flags = int_ptr(n);
+                    }
+                } else if(prop_type == Prop_Real) {
+                    auto float_ptr = ps->getAsFloatProperty(prop);
+                    float_ptr(n) = std::stod(in0);
+                } else {
+                    std::cerr << "read_particle_data(): Invalid property type!" << std::endl;
+                    return 0;
+                }
+
+                i++;
+            }
+
+            if(within_domain || flags & (FLAGS_INFINITE | FLAGS_FIXED | FLAGS_GLOBAL)) {
+                shape_ptr(n++) = shape_id;
+            }
+        }
+
+        in_file.close();
+    }
+
+    return n;
+}
+
+/*
+size_t read_feature_data(PairsRuntime *ps, const char *filename, const int feature_id, const property_t properties[], size_t nprops) {
+    std::ifstream in_file(filename, std::ifstream::in);
+    std::string line;
+
+    if(in_file.is_open()) {
+        while(std::getline(in_file, line)) {
+            std::stringstream line_stream(line);
+            std::string istr, jstr, in0;
+            std::getline(line_stream, istr, ',');
+            std::getline(line_stream, jstr, ',');
+            int i = std::stoi(istr);
+            int j = std::stoi(jstr);
+
+            while(std::getline(line_stream, in0, ',')) {
+                property_t p_id = properties[i];
+                auto prop = ps->getProperty(p_id);
+                auto prop_type = prop.getType();
+
+                if(prop_type == Prop_Vector) {
+                    auto vector_ptr = ps->getAsVectorFeatureProperty(prop);
+                    std::string in1, in2;
+                    std::getline(line_stream, in1, ',');
+                    std::getline(line_stream, in2, ',');
+                    real_t x = std::stod(in0);
+                    real_t y = std::stod(in1);
+                    real_t z = std::stod(in2);
+                    vector_ptr(i, j, 0) = x;
+                    vector_ptr(i, j, 1) = y;
+                    vector_ptr(i, j, 2) = z;
+                } else if(prop_type == Prop_Integer) {
+                    auto int_ptr = ps->getAsIntegerFeatureProperty(prop);
+                    int_ptr(i, j) = std::stoi(in0);
+                } else if(prop_type == Prop_Real) {
+                    auto float_ptr = ps->getAsFloatFeatureProperty(prop);
+                    float_ptr(i, j) = std::stod(in0);
+                } else {
+                    std::cerr << "read_feature_data(): Invalid property type!" << std::endl;
+                    return 0;
+                }
+            }
+        }
+
+        in_file.close();
+    }
+
+    return n;
+}
+*/
+
+}
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index cae3362b0fce5997c516fbd1be9220953ca4d16b..ca35f316fc2d6a6044f5317bca9cb222632abb84 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -1,8 +1,3 @@
-#include <iostream>
-#include <string.h>
-#include <fstream>
-#include <sstream>
-//---
 #include "pairs.hpp"
 #include "pairs_common.hpp"
 
@@ -10,158 +5,10 @@
 
 namespace pairs {
 
-void read_grid_data(PairsRuntime *ps, const char *filename, real_t *grid_buffer) {
-    std::ifstream in_file(filename, std::ifstream::in);
-    std::string line;
+void read_grid_data(PairsRuntime *ps, const char *filename, real_t *grid_buffer);
 
-    if(in_file.is_open()) {
-        std::getline(in_file, line);
-        std::stringstream line_stream(line);
-        std::string in0;
-        int i = 0;
-
-        while(std::getline(line_stream, in0, ',')) {
-            //PAIRS_ASSERT(i < ndims * 2);
-            grid_buffer[i] = std::stod(in0);
-            i++;
-        }
-
-        in_file.close();
-    }
-}
-
-size_t read_particle_data(PairsRuntime *ps, const char *filename, const property_t properties[], size_t nprops, int shape_id, int start) {
-    std::ifstream in_file(filename, std::ifstream::in);
-    std::string line;
-    auto shape_ptr = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
-    size_t n = start;
-
-    if(in_file.is_open()) {
-        //std::getline(in_file, line);
-        while(std::getline(in_file, line)) {
-            std::stringstream line_stream(line);
-            std::string in0;
-            int within_domain = 1;
-            int i = 0;
-            int flags = 0;
-
-            while(std::getline(line_stream, in0, ',')) {
-                property_t p_id = properties[i];
-                auto prop = ps->getProperty(p_id);
-                auto prop_type = prop.getType();
-
-                if(prop_type == Prop_Vector) {
-                    auto vector_ptr = ps->getAsVectorProperty(prop);
-                    std::string in1, in2;
-                    std::getline(line_stream, in1, ',');
-                    std::getline(line_stream, in2, ',');
-                    real_t x = std::stod(in0);
-                    real_t y = std::stod(in1);
-                    real_t z = std::stod(in2);
-                    vector_ptr(n, 0) = x;
-                    vector_ptr(n, 1) = y;
-                    vector_ptr(n, 2) = z;
-
-                    if(prop.getName() == "position") {
-                        within_domain = ps->getDomainPartitioner()->isWithinSubdomain(x, y, z);
-                    }
-                } else if(prop_type == Prop_Matrix) {
-                    auto matrix_ptr = ps->getAsMatrixProperty(prop);
-                    constexpr int nelems = 9;
-                    std::string in_buf;
-
-                    matrix_ptr(n, 0) = std::stod(in0);
-                    for(int i = 1; i < nelems; i++) {
-                        std::getline(line_stream, in_buf, ',');
-                        matrix_ptr(n, i) = std::stod(in_buf);
-                    }
-                } else if(prop_type == Prop_Quaternion) {
-                    auto quat_ptr = ps->getAsQuaternionProperty(prop);
-                    constexpr int nelems = 4;
-                    std::string in_buf;
-
-                    quat_ptr(n, 0) = std::stod(in0);
-                    for(int i = 1; i < nelems; i++) {
-                        std::getline(line_stream, in_buf, ',');
-                        quat_ptr(n, i) = std::stod(in_buf);
-                    }
-                } else if(prop_type == Prop_Integer) {
-                    auto int_ptr = ps->getAsIntegerProperty(prop);
-                    int_ptr(n) = std::stoi(in0);
-
-                    if(prop.getName() == "flags") {
-                        flags = int_ptr(n);
-                    }
-                } else if(prop_type == Prop_Real) {
-                    auto float_ptr = ps->getAsFloatProperty(prop);
-                    float_ptr(n) = std::stod(in0);
-                } else {
-                    std::cerr << "read_particle_data(): Invalid property type!" << std::endl;
-                    return 0;
-                }
-
-                i++;
-            }
-
-            if(within_domain || flags & (FLAGS_INFINITE | FLAGS_FIXED | FLAGS_GLOBAL)) {
-                shape_ptr(n++) = shape_id;
-            }
-        }
-
-        in_file.close();
-    }
-
-    return n;
-}
-
-/*
-size_t read_feature_data(PairsRuntime *ps, const char *filename, const int feature_id, const property_t properties[], size_t nprops) {
-    std::ifstream in_file(filename, std::ifstream::in);
-    std::string line;
-
-    if(in_file.is_open()) {
-        while(std::getline(in_file, line)) {
-            std::stringstream line_stream(line);
-            std::string istr, jstr, in0;
-            std::getline(line_stream, istr, ',');
-            std::getline(line_stream, jstr, ',');
-            int i = std::stoi(istr);
-            int j = std::stoi(jstr);
-
-            while(std::getline(line_stream, in0, ',')) {
-                property_t p_id = properties[i];
-                auto prop = ps->getProperty(p_id);
-                auto prop_type = prop.getType();
-
-                if(prop_type == Prop_Vector) {
-                    auto vector_ptr = ps->getAsVectorFeatureProperty(prop);
-                    std::string in1, in2;
-                    std::getline(line_stream, in1, ',');
-                    std::getline(line_stream, in2, ',');
-                    real_t x = std::stod(in0);
-                    real_t y = std::stod(in1);
-                    real_t z = std::stod(in2);
-                    vector_ptr(i, j, 0) = x;
-                    vector_ptr(i, j, 1) = y;
-                    vector_ptr(i, j, 2) = z;
-                } else if(prop_type == Prop_Integer) {
-                    auto int_ptr = ps->getAsIntegerFeatureProperty(prop);
-                    int_ptr(i, j) = std::stoi(in0);
-                } else if(prop_type == Prop_Real) {
-                    auto float_ptr = ps->getAsFloatFeatureProperty(prop);
-                    float_ptr(i, j) = std::stod(in0);
-                } else {
-                    std::cerr << "read_feature_data(): Invalid property type!" << std::endl;
-                    return 0;
-                }
-            }
-        }
-
-        in_file.close();
-    }
-
-    return n;
-}
-*/
+size_t read_particle_data(
+    PairsRuntime *ps, const char *filename, const property_t properties[],
+    size_t nprops, int shape_id, int start);
 
 }
diff --git a/runtime/stats.cpp b/runtime/stats.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d66c34162937c36d2827be8d1bc8c3852461013c
--- /dev/null
+++ b/runtime/stats.cpp
@@ -0,0 +1,29 @@
+#include "pairs.hpp"
+
+namespace pairs {
+
+void print_stats(PairsRuntime *ps, int nlocal, int nghost) {
+    int min_nlocal = nlocal;
+    int max_nlocal = nlocal;
+    int min_nghost = nghost;
+    int max_nghost = nghost;
+    int nglobal;
+
+    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+        MPI_Allreduce(&nlocal, &nglobal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+        min_nlocal = nglobal;
+        MPI_Allreduce(&nlocal, &nglobal, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+        max_nlocal = nglobal;
+        MPI_Allreduce(&nghost, &nglobal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+        min_nghost = nglobal;
+        MPI_Allreduce(&nghost, &nglobal, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+        max_nghost = nglobal;
+    }
+
+    if(ps->getDomainPartitioner()->getRank() == 0) {
+        std::cout << "Number of local particles: " << min_nlocal << " / " << max_nlocal << std::endl;
+        std::cout << "Number of ghost particles: " << min_nghost << " / " << max_nghost << std::endl;
+    }
+}
+
+}
diff --git a/runtime/stats.hpp b/runtime/stats.hpp
index 0bb9e9bb8c606f97a57bba0af51b0d82ee37277d..e6c51c306f83b347534963aac7c24f7a42f1b58d 100644
--- a/runtime/stats.hpp
+++ b/runtime/stats.hpp
@@ -2,32 +2,8 @@
 
 #pragma once
 
-using namespace std;
-
 namespace pairs {
 
-void print_stats(PairsRuntime *ps, int nlocal, int nghost) {
-    int min_nlocal = nlocal;
-    int max_nlocal = nlocal;
-    int min_nghost = nghost;
-    int max_nghost = nghost;
-    int nglobal;
-
-    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
-        MPI_Allreduce(&nlocal, &nglobal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
-        min_nlocal = nglobal;
-        MPI_Allreduce(&nlocal, &nglobal, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
-        max_nlocal = nglobal;
-        MPI_Allreduce(&nghost, &nglobal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
-        min_nghost = nglobal;
-        MPI_Allreduce(&nghost, &nglobal, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
-        max_nghost = nglobal;
-    }
-
-    if(ps->getDomainPartitioner()->getRank() == 0) {
-        std::cout << "Number of local particles: " << min_nlocal << " / " << max_nlocal << std::endl;
-        std::cout << "Number of ghost particles: " << min_nghost << " / " << max_nghost << std::endl;
-    }
-}
+void print_stats(PairsRuntime *ps, int nlocal, int nghost);
 
 }
diff --git a/runtime/thermo.cpp b/runtime/thermo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..044f6d3639f8db286099a19b112b67b065f2cfa0
--- /dev/null
+++ b/runtime/thermo.cpp
@@ -0,0 +1,101 @@
+#include <iostream>
+#include <math.h>
+#include <mpi.h>
+//---
+#include "pairs.hpp"
+
+namespace pairs {
+
+double compute_thermo(
+    PairsRuntime *ps, int nlocal, double xprd, double yprd, double zprd, int print) {
+
+    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
+    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+    int natoms = nlocal;
+
+    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+        int global_natoms;
+        MPI_Allreduce(&natoms, &global_natoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+        natoms = global_natoms;
+    }
+
+    const double mvv2e = 1.0;
+    const double dof_boltz = (natoms * 3 - 3);
+    const double t_scale = mvv2e / dof_boltz;
+    const double p_scale = 1.0 / 3 / xprd / yprd / zprd;
+    //const double e_scale = 0.5;
+    double t = 0.0, p;
+
+    ps->copyPropertyToHost(masses, ReadOnly);
+    ps->copyPropertyToHost(velocities, ReadOnly);
+
+    for(int i = 0; i < nlocal; i++) {
+        t += masses(i) * (  velocities(i, 0) * velocities(i, 0) +
+                            velocities(i, 1) * velocities(i, 1) +
+                            velocities(i, 2) * velocities(i, 2)   );
+    }
+
+    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+        double global_t;
+        MPI_Allreduce(&t, &global_t, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        t = global_t;
+    }
+
+    t = t * t_scale;
+    if(print == 1 && ps->getDomainPartitioner()->getRank() == 0) {
+        p = (t * dof_boltz) * p_scale;
+        std::cout << t << "\t" << p << std::endl;
+    }
+
+    return t;
+}
+
+void adjust_thermo(
+    PairsRuntime *ps, int nlocal, double xprd, double yprd, double zprd, double temp) {
+
+    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+    double vxtot = 0.0;
+    double vytot = 0.0;
+    double vztot = 0.0;
+    double tmp;
+    int natoms = nlocal;
+
+    for(int i = 0; i < nlocal; i++) {
+        vxtot += velocities(i, 0);
+        vytot += velocities(i, 1);
+        vztot += velocities(i, 2);
+    }
+
+    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+        int global_natoms;
+        MPI_Allreduce(&natoms, &global_natoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+        natoms = global_natoms;
+        MPI_Allreduce(&vxtot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        vxtot = tmp / natoms;
+        MPI_Allreduce(&vytot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        vytot = tmp / natoms;
+        MPI_Allreduce(&vztot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        vztot = tmp / natoms;
+    } else {
+        vxtot /= natoms;
+        vytot /= natoms;
+        vztot /= natoms;
+    }
+
+    for(int i = 0; i < nlocal; i++) {
+        velocities(i, 0) -= vxtot;
+        velocities(i, 1) -= vytot;
+        velocities(i, 2) -= vztot;
+    }
+
+    double t = pairs::compute_thermo(ps, nlocal, xprd, yprd, zprd, 0);
+    double factor = sqrt(temp / t);
+
+    for(int i = 0; i < nlocal; i++) {
+        velocities(i, 0) *= factor;
+        velocities(i, 1) *= factor;
+        velocities(i, 2) *= factor;
+    }
+}
+
+}
diff --git a/runtime/thermo.hpp b/runtime/thermo.hpp
index 71385924b632f01acaf532c94aaf13dea3cff909..6902b007603478ec2f4ca23136c58a88b025eef1 100644
--- a/runtime/thermo.hpp
+++ b/runtime/thermo.hpp
@@ -1,99 +1,13 @@
-#include <iostream>
-#include <math.h>
-#include <mpi.h>
-//---
 #include "pairs.hpp"
 
 #pragma once
 
 namespace pairs {
 
-double compute_thermo(PairsRuntime *ps, int nlocal, double xprd, double yprd, double zprd, int print) {
-    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
-    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
-    int natoms = nlocal;
+double compute_thermo(
+    PairsRuntime *ps, int nlocal, double xprd, double yprd, double zprd, int print);
 
-    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
-        int global_natoms;
-        MPI_Allreduce(&natoms, &global_natoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-        natoms = global_natoms;
-    }
-
-    const double mvv2e = 1.0;
-    const double dof_boltz = (natoms * 3 - 3);
-    const double t_scale = mvv2e / dof_boltz;
-    const double p_scale = 1.0 / 3 / xprd / yprd / zprd;
-    //const double e_scale = 0.5;
-    double t = 0.0, p;
-
-    ps->copyPropertyToHost(masses, ReadOnly);
-    ps->copyPropertyToHost(velocities, ReadOnly);
-
-    for(int i = 0; i < nlocal; i++) {
-        t += masses(i) * (  velocities(i, 0) * velocities(i, 0) +
-                            velocities(i, 1) * velocities(i, 1) +
-                            velocities(i, 2) * velocities(i, 2)   );
-    }
-
-    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
-        double global_t;
-        MPI_Allreduce(&t, &global_t, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-        t = global_t;
-    }
-
-    t = t * t_scale;
-    if(print == 1 && ps->getDomainPartitioner()->getRank() == 0) {
-        p = (t * dof_boltz) * p_scale;
-        std::cout << t << "\t" << p << std::endl;
-    }
-
-    return t;
-}
-
-void adjust_thermo(PairsRuntime *ps, int nlocal, double xprd, double yprd, double zprd, double temp) {
-    auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
-    double vxtot = 0.0;
-    double vytot = 0.0;
-    double vztot = 0.0;
-    double tmp;
-    int natoms = nlocal;
-
-    for(int i = 0; i < nlocal; i++) {
-        vxtot += velocities(i, 0);
-        vytot += velocities(i, 1);
-        vztot += velocities(i, 2);
-    }
-
-    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
-        int global_natoms;
-        MPI_Allreduce(&natoms, &global_natoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-        natoms = global_natoms;
-        MPI_Allreduce(&vxtot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-        vxtot = tmp / natoms;
-        MPI_Allreduce(&vytot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-        vytot = tmp / natoms;
-        MPI_Allreduce(&vztot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-        vztot = tmp / natoms;
-    } else {
-        vxtot /= natoms;
-        vytot /= natoms;
-        vztot /= natoms;
-    }
-
-    for(int i = 0; i < nlocal; i++) {
-        velocities(i, 0) -= vxtot;
-        velocities(i, 1) -= vytot;
-        velocities(i, 2) -= vztot;
-    }
-
-    double t = pairs::compute_thermo(ps, nlocal, xprd, yprd, zprd, 0);
-    double factor = sqrt(temp / t);
-
-    for(int i = 0; i < nlocal; i++) {
-        velocities(i, 0) *= factor;
-        velocities(i, 1) *= factor;
-        velocities(i, 2) *= factor;
-    }
-}
+void adjust_thermo(
+    PairsRuntime *ps, int nlocal, double xprd, double yprd, double zprd, double temp);
 
 }
diff --git a/runtime/timing.cpp b/runtime/timing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0068d8117b89d2529e3f8dc10b24e1d2483672e8
--- /dev/null
+++ b/runtime/timing.cpp
@@ -0,0 +1,23 @@
+#include "pairs.hpp"
+
+using namespace std;
+
+namespace pairs {
+
+void register_timer(PairsRuntime *ps, int id, std::string name) {
+    ps->getTimers()->add(id, name);
+}
+
+void start_timer(PairsRuntime *ps, int id) {
+    ps->getTimers()->start(id);
+}
+
+void stop_timer(PairsRuntime *ps, int id) {
+    ps->getTimers()->stop(id);
+}
+
+void print_timers(PairsRuntime *ps) {
+    ps->printTimers();
+}
+
+}
diff --git a/runtime/timing.hpp b/runtime/timing.hpp
index 77a799dbf20266570d4e48245eb203e7dca51db7..f7544603549232cce89c00e7071948da32511693 100644
--- a/runtime/timing.hpp
+++ b/runtime/timing.hpp
@@ -6,20 +6,9 @@ using namespace std;
 
 namespace pairs {
 
-void register_timer(PairsRuntime *ps, int id, std::string name) {
-    ps->getTimers()->add(id, name);
-}
-
-void start_timer(PairsRuntime *ps, int id) {
-    ps->getTimers()->start(id);
-}
-
-void stop_timer(PairsRuntime *ps, int id) {
-    ps->getTimers()->stop(id);
-}
-
-void print_timers(PairsRuntime *ps) {
-    ps->printTimers();
-}
+void register_timer(PairsRuntime *ps, int id, std::string name);
+void start_timer(PairsRuntime *ps, int id);
+void stop_timer(PairsRuntime *ps, int id);
+void print_timers(PairsRuntime *ps);
 
 }
diff --git a/runtime/vtk.cpp b/runtime/vtk.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cc47050082308d6171ba57df652dfdf40408af78
--- /dev/null
+++ b/runtime/vtk.cpp
@@ -0,0 +1,88 @@
+#include <iomanip>
+#include <iostream>
+#include <fstream>
+//---
+#include "pairs.hpp"
+
+namespace pairs {
+
+void vtk_write_data(
+    PairsRuntime *ps, const char *filename, int start, int end, int timestep, int frequency) {
+
+    std::string output_filename(filename);
+    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
+    auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
+    auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
+    const int prec = 8;
+    int n = end - start;
+    std::ostringstream filename_oss;
+
+    if(frequency != 0 && timestep % frequency != 0) {
+        return;
+    }
+
+    filename_oss << filename << "_";
+    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+        filename_oss << "r" << ps->getDomainPartitioner()->getRank() << "_";
+    }
+
+    filename_oss << timestep << ".vtk";
+    std::ofstream out_file(filename_oss.str());
+
+    ps->copyPropertyToHost(masses, ReadOnly);
+    ps->copyPropertyToHost(positions, ReadOnly);
+    ps->copyPropertyToHost(flags, ReadOnly);
+
+    for(int i = start; i < end; i++) {
+        if(flags(i) & FLAGS_INFINITE) {
+            n--;
+        }
+    }
+
+    if(out_file.is_open()) {
+        out_file << "# vtk DataFile Version 2.0\n";
+        out_file << "Particle data\n";
+        out_file << "ASCII\n";
+        out_file << "DATASET UNSTRUCTURED_GRID\n";
+        out_file << "POINTS " << n << " double\n";
+
+        for(int i = start; i < end; i++) {
+            if(!(flags(i) & FLAGS_INFINITE)) {
+                out_file << std::fixed << std::setprecision(prec) << positions(i, 0) << " ";
+                out_file << std::fixed << std::setprecision(prec) << positions(i, 1) << " ";
+                out_file << std::fixed << std::setprecision(prec) << positions(i, 2) << "\n";
+            }
+        }
+
+        out_file << "\n\n";
+        out_file << "CELLS " << n << " " << (n * 2) << "\n";
+        for(int i = start; i < end; i++) {
+            if(!(flags(i) & FLAGS_INFINITE)) {
+                out_file << "1 " << (i - start) << "\n";
+            }
+        }
+
+        out_file << "\n\n";
+        out_file << "CELL_TYPES " << n << "\n";
+        for(int i = start; i < end; i++) {
+            if(!(flags(i) & FLAGS_INFINITE)) {
+                out_file << "1\n";
+            }
+        }
+
+        out_file << "\n\n";
+        out_file << "POINT_DATA " << n << "\n";
+        out_file << "SCALARS mass double\n";
+        out_file << "LOOKUP_TABLE default\n";
+        for(int i = start; i < end; i++) {
+            if(!(flags(i) & FLAGS_INFINITE)) {
+                out_file << std::fixed << std::setprecision(prec) << masses(i) << "\n";
+            }
+        }
+
+        out_file << "\n\n";
+        out_file.close();
+    }
+}
+
+}
diff --git a/runtime/vtk.hpp b/runtime/vtk.hpp
index 181c2181f3a34a2f71b72545bef22019d65a5a1b..2bc766ca3f065bbf3f96b528de33c8b284e8bf87 100644
--- a/runtime/vtk.hpp
+++ b/runtime/vtk.hpp
@@ -1,88 +1,10 @@
-#include <iomanip>
-#include <iostream>
-#include <fstream>
-//---
 #include "pairs.hpp"
 
 #pragma once
 
 namespace pairs {
 
-void vtk_write_data(PairsRuntime *ps, const char *filename, int start, int end, int timestep, int frequency) {
-    std::string output_filename(filename);
-    auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
-    auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
-    auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
-    const int prec = 8;
-    int n = end - start;
-    std::ostringstream filename_oss;
-
-    if(frequency != 0 && timestep % frequency != 0) {
-        return;
-    }
-
-    filename_oss << filename << "_";
-    if(ps->getDomainPartitioner()->getWorldSize() > 1) {
-        filename_oss << "r" << ps->getDomainPartitioner()->getRank() << "_";
-    }
-
-    filename_oss << timestep << ".vtk";
-    std::ofstream out_file(filename_oss.str());
-
-    ps->copyPropertyToHost(masses, ReadOnly);
-    ps->copyPropertyToHost(positions, ReadOnly);
-    ps->copyPropertyToHost(flags, ReadOnly);
-
-    for(int i = start; i < end; i++) {
-        if(flags(i) & FLAGS_INFINITE) {
-            n--;
-        }
-    }
-
-    if(out_file.is_open()) {
-        out_file << "# vtk DataFile Version 2.0\n";
-        out_file << "Particle data\n";
-        out_file << "ASCII\n";
-        out_file << "DATASET UNSTRUCTURED_GRID\n";
-        out_file << "POINTS " << n << " double\n";
-
-        for(int i = start; i < end; i++) {
-            if(!(flags(i) & FLAGS_INFINITE)) {
-                out_file << std::fixed << std::setprecision(prec) << positions(i, 0) << " ";
-                out_file << std::fixed << std::setprecision(prec) << positions(i, 1) << " ";
-                out_file << std::fixed << std::setprecision(prec) << positions(i, 2) << "\n";
-            }
-        }
-
-        out_file << "\n\n";
-        out_file << "CELLS " << n << " " << (n * 2) << "\n";
-        for(int i = start; i < end; i++) {
-            if(!(flags(i) & FLAGS_INFINITE)) {
-                out_file << "1 " << (i - start) << "\n";
-            }
-        }
-
-        out_file << "\n\n";
-        out_file << "CELL_TYPES " << n << "\n";
-        for(int i = start; i < end; i++) {
-            if(!(flags(i) & FLAGS_INFINITE)) {
-                out_file << "1\n";
-            }
-        }
-
-        out_file << "\n\n";
-        out_file << "POINT_DATA " << n << "\n";
-        out_file << "SCALARS mass double\n";
-        out_file << "LOOKUP_TABLE default\n";
-        for(int i = start; i < end; i++) {
-            if(!(flags(i) & FLAGS_INFINITE)) {
-                out_file << std::fixed << std::setprecision(prec) << masses(i) << "\n";
-            }
-        }
-
-        out_file << "\n\n";
-        out_file.close();
-    }
-}
+void vtk_write_data(
+    PairsRuntime *ps, const char *filename, int start, int end, int timestep, int frequency);
 
 }