             cudaDeviceSynchronize();
 		}
 
-        // Utility function for the SIB Optimizer
-        __global__ void cu_get_random_vectorfield(scalar epsilon, Vector3 * xi, size_t N)
+        void get_random_vector(std::uniform_real_distribution<scalar> & distribution, std::mt19937 & prng, Vector3 & vec)
+        {
+            for (int dim = 0; dim < 3; ++dim)
+            {
+                vec[dim] = distribution(prng);
+            }
+        }
+        
+        __global__ void cu_get_random_vectorfield(Vector3 * xi, size_t N)
         {
             unsigned long long subsequence = 0;
             unsigned long long offset= 0;

                curand_init(idx,subsequence,offset,&state);
                for (int dim=0;dim<3; ++dim)
                {
                    xi[idx][dim] = epsilon*(llroundf(curand_uniform(&state))*2-1);
                }
            }
        }
		void get_random_vectorfield(const Data::Spin_System & sys, scalar epsilon, vectorfield & xi)
		{
            int n = xi.size();
            cu_get_random_vectorfield<<<(n+1023)/1024, 1024>>>(epsilon, xi.data(), n);
            cudaDeviceSynchronize();
		}

        int neigh_cu_get_pair_j(const int * boundary_conditions, const int * n_cells, int N, int ispin, Neighbour neigh)
        {

            cudaDeviceSynchronize();
        }

        __global__ void cu_add_c_a2_mask(scalar c, const Vector3 * a, Vector3 * out, const int * mask, size_t N)
        {
            int idx = blockIdx.x * blockDim.x + threadIdx.x;
            if(idx < N)
            {
                out[idx] += c*mask[idx]*a[idx];
            }
        }
        // out[i] += c*a[i]
        void add_c_a(const scalar & c, const vectorfield & a, vectorfield & out, const intfield & mask)
        {
            int n = out.size();
            cu_add_c_a2_mask<<<(n+1023)/1024, 1024>>>(c, a.data(), out.data(), mask.data(), n);
            cudaDeviceSynchronize();
        }

        __global__ void cu_add_c_a3(const scalar * c, const Vector3 * a, Vector3 * out, size_t N)
        {
            int idx = blockIdx.x * blockDim.x + threadIdx.x;
