diff --git a/examples/benchmarks/sd_static.cpp b/examples/benchmarks/sd_static.cpp
index aad480137a1ece7601baf1e170fdcbe4a7bd0188..96bed6381b398063f815dccf3d617184280c4c90 100644
--- a/examples/benchmarks/sd_static.cpp
+++ b/examples/benchmarks/sd_static.cpp
@@ -4,6 +4,21 @@
 
 #include "spring_dashpot.hpp"
 
+void randomaize_indices(PairsAccessor ac){
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dist(0, ac.nlocal() - 1);
+
+    ac.syncPosition(PairsAccessor::Host);
+    for (int i=0; i<ac.nlocal(); ++i){
+        int j = dist(gen);
+        auto tmp = ac.getPosition(i);
+        ac.setPosition(i, ac.getPosition(j));
+        ac.setPosition(j, tmp);
+    }
+    ac.syncPosition(PairsAccessor::Host);
+}
+
 int main(int argc, char **argv) {
     if(argc!=5){
         std::cerr << "4 args are required: Domain size (i.e. number of particles) in x,y,z and #timesteps." << std::endl;
@@ -15,6 +30,8 @@ int main(int argc, char **argv) {
 
     auto pairs_sim = std::make_shared<PairsSimulation>();
     pairs_sim->initialize();
+    
+    PairsAccessor ac(pairs_sim.get());
 
     auto pairs_runtime = pairs_sim->getPairsRuntime();
 
@@ -44,6 +61,9 @@ int main(int argc, char **argv) {
     pairs_sim->setInteractionRadius(cell_width);
     pairs_sim->updateDomain();
     
+    randomaize_indices(ac);
+    pairs_sim->reneighbor();
+
     // Inertia update is required for euler updates to be valid (but particles remain stationary)
     double dt = 0.001;  // Arbitrary
     
@@ -89,5 +109,6 @@ int main(int argc, char **argv) {
     // pairs::vtk_write_data(pairs_runtime, "output/local_spheres", 0, pairs_sim->nlocal(), 0);
     // pairs::vtk_write_data(pairs_runtime, "output/ghost_spheres", pairs_sim->nlocal(), pairs_sim->size(), 0);
     
+    ac.end();
     pairs_sim->end();
 }
\ No newline at end of file
diff --git a/examples/benchmarks/sd_static_triangular.cpp b/examples/benchmarks/sd_static_triangular.cpp
index 284b8aa8c6777ad0cd64b0c2fd7658ad37077284..7a05fe564f890c2e06a27f77eece8599575e10ed 100644
--- a/examples/benchmarks/sd_static_triangular.cpp
+++ b/examples/benchmarks/sd_static_triangular.cpp
@@ -4,6 +4,21 @@
 
 #include "spring_dashpot_no_pbc.hpp"
 
+void randomaize_indices(PairsAccessor ac){
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dist(0, ac.nlocal() - 1);
+
+    ac.syncPosition(PairsAccessor::Host);
+    for (int i=0; i<ac.nlocal(); ++i){
+        int j = dist(gen);
+        auto tmp = ac.getPosition(i);
+        ac.setPosition(i, ac.getPosition(j));
+        ac.setPosition(j, tmp);
+    }
+    ac.syncPosition(PairsAccessor::Host);
+}
+
 template<typename Type>
 void print_global_stats(std::string name, Type value, MPI_Datatype mpi_type, MPI_Comm comm){
     int rank, world_size;
@@ -60,8 +75,10 @@ int main(int argc, char **argv) {
     auto pairs_sim = std::make_shared<PairsSimulation>();
     pairs_sim->initialize();
 
+    PairsAccessor ac(pairs_sim.get());
+    
     auto pairs_runtime = pairs_sim->getPairsRuntime();
-
+    
     bool load_balanced = (std::stoi(argv[5]) != 0);
     pairs_runtime->initDomain(&argc, &argv, 0, 0, 0, domain_size[0], domain_size[1], domain_size[2], load_balanced);
 
@@ -110,6 +127,9 @@ int main(int argc, char **argv) {
     // Rebalance
     pairs_sim->updateDomain();
 
+    randomaize_indices(ac);
+    pairs_sim->reneighbor();
+
     // Stats
     // ------------------------------------------------------------------------------
     int rank = pairs_sim->rank();
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 31540997c862f9458ff7e3e59dae7eeaa6c853d6..62ecdd8599a07c5e3c7d10296b430ca80021b6bc 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -2,6 +2,7 @@
 #include <cmath>
 #include <memory>
 #include <vector>
+#include <array>
 //---
 #include "allocate.hpp"
 #include "array.hpp"