diff --git a/core/include/engine/Indexing.hpp b/core/include/engine/Indexing.hpp
new file mode 100644
index 000000000..cbff4f17a
--- /dev/null
+++ b/core/include/engine/Indexing.hpp
@@ -0,0 +1,450 @@
+#pragma once
+#ifndef SPIRIT_CORE_ENGINE_INDEXING_HPP
+#define SPIRIT_CORE_ENGINE_INDEXING_HPP
+
+#include <engine/Vectormath_Defines.hpp>
+
+namespace Engine
+{
+namespace Vectormath
+{
+
+// Note: translations must lie within bounds of n_cells
+inline int
+idx_from_translations( const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations )
+{
+    auto & Na = n_cells[0];
+    auto & Nb = n_cells[1];
+    auto & Nc = n_cells[2];
+    auto & N  = n_cell_atoms;
+
+    auto & da = translations[0];
+    auto & db = translations[1];
+    auto & dc = translations[2];
+
+    return da * N + db * N * Na + dc * N * Na * Nb;
+}
+
+#ifndef SPIRIT_USE_CUDA
+
+// Get the linear index in a n-D array where tupel contains the components in n-dimensions from fatest to slowest
+// varying and maxVal is the extent in every dimension
+inline int idx_from_tupel( const field<int> & tupel, const field<int> & maxVal )
+{
+    int idx  = 0;
+    int mult = 1;
+    for( int i = 0; i < tupel.size(); i++ )
+    {
+        idx += mult * tupel[i];
+        mult *= maxVal[i];
+    }
+    return idx;
+}
+
+// reverse of idx_from_tupel
+inline void tupel_from_idx( int & idx, field<int> & tupel, field<int> & maxVal )
+{
+    int idx_diff = idx;
+    int div      = 1;
+    for( int i = 0; i < tupel.size() - 1; i++ )
+        div *= maxVal[i];
+    for( int i = tupel.size() - 1; i > 0; i-- )
+    {
+        tupel[i] = idx_diff / div;
+        idx_diff -= tupel[i] * div;
+        div /= maxVal[i - 1];
+    }
+    tupel[0] = idx_diff / div;
+}
+
+inline int idx_from_translations(
+    const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations_i,
+    const std::array<int, 3> & translations )
+{
+    auto & Na = n_cells[0];
+    auto & Nb = n_cells[1];
+    auto & Nc = n_cells[2];
+    auto & N  = n_cell_atoms;
+
+    int da = translations_i[0] + translations[0];
+    int db = translations_i[1] + translations[1];
+    int dc = translations_i[2] + translations[2];
+
+    if( translations[0] < 0 )
+        da += N * Na;
+    if( translations[1] < 0 )
+        db += N * Na * Nb;
+    if( translations[2] < 0 )
+        dc += N * Na * Nb * Nc;
+
+    return ( da % Na ) * N + ( db % Nb ) * N * Na + ( dc % Nc ) * N * Na * Nb;
+}
+
+inline bool boundary_conditions_fulfilled(
+    const intfield & n_cells, const intfield & boundary_conditions, const std::array<int, 3> & translations_i,
+    const std::array<int, 3> & translations_j )
+{
+    int da = translations_i[0] + translations_j[0];
+    int db = translations_i[1] + translations_j[1];
+    int dc = translations_i[2] + translations_j[2];
+    return (
+        ( boundary_conditions[0] || ( 0 <= da && da < n_cells[0] ) )
+        && ( boundary_conditions[1] || ( 0 <= db && db < n_cells[1] ) )
+        && ( boundary_conditions[2] || ( 0 <= dc && dc < n_cells[2] ) ) );
+}
+
+#endif
+#ifdef SPIRIT_USE_CUDA
+
+// Get the linear index in a n-D array where tupel contains the components in n-dimensions from fatest to slowest
+// varying and maxVal is the extent in every dimension
+inline __device__ int cu_idx_from_tupel( field<int> & tupel, field<int> & maxVal )
+{
+    int idx  = 0;
+    int mult = 1;
+    for( int i = 0; i < tupel.size(); i++ )
+    {
+        idx += mult * tupel[i];
+        mult *= maxVal[i];
+    }
+    return idx;
+}
+
+// reverse of idx_from_tupel
+inline __device__ void cu_tupel_from_idx( int & idx, int * tupel, int * maxVal, int n )
+{
+    int idx_diff = idx;
+    int div      = 1;
+    for( int i = 0; i < n - 1; i++ )
+        div *= maxVal[i];
+    for( int i = n - 1; i > 0; i-- )
+    {
+        tupel[i] = idx_diff / div;
+        idx_diff -= tupel[i] * div;
+        div /= maxVal[i - 1];
+    }
+    tupel[0] = idx_diff / div;
+}
+
+inline void tupel_from_idx( int & idx, field<int> & tupel, const field<int> & maxVal )
+{
+    int idx_diff = idx;
+    int div      = 1;
+    for( int i = 0; i < maxVal.size() - 1; i++ )
+        div *= maxVal[i];
+    for( int i = maxVal.size() - 1; i > 0; i-- )
+    {
+        tupel[i] = idx_diff / div;
+        idx_diff -= tupel[i] * div;
+        div /= maxVal[i - 1];
+    }
+    tupel[0] = idx_diff / div;
+}
+
+inline int idx_from_translations(
+    const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations_i,
+    const int translations[3] )
+{
+    int Na = n_cells[0];
+    int Nb = n_cells[1];
+    int Nc = n_cells[2];
+    int N  = n_cell_atoms;
+
+    int da = translations_i[0] + translations[0];
+    int db = translations_i[1] + translations[1];
+    int dc = translations_i[2] + translations[2];
+
+    if( translations[0] < 0 )
+        da += N * Na;
+    if( translations[1] < 0 )
+        db += N * Na * Nb;
+    if( translations[2] < 0 )
+        dc += N * Na * Nb * Nc;
+
+    int idx = ( da % Na ) * N + ( db % Nb ) * N * Na + ( dc % Nc ) * N * Na * Nb;
+
+    return idx;
+}
+
+inline bool boundary_conditions_fulfilled(
+    const intfield & n_cells, const intfield & boundary_conditions, const std::array<int, 3> & translations_i,
+    const int translations_j[3] )
+{
+    int da = translations_i[0] + translations_j[0];
+    int db = translations_i[1] + translations_j[1];
+    int dc = translations_i[2] + translations_j[2];
+    return (
+        ( boundary_conditions[0] || ( 0 <= da && da < n_cells[0] ) )
+        && ( boundary_conditions[1] || ( 0 <= db && db < n_cells[1] ) )
+        && ( boundary_conditions[2] || ( 0 <= dc && dc < n_cells[2] ) ) );
+}
+
+__inline__ __device__ bool cu_check_atom_type( int atom_type )
+{
+#ifdef SPIRIT_ENABLE_DEFECTS
+    // If defects are enabled we check for
+    //		vacancies (type < 0)
+    if( atom_type >= 0 )
+        return true;
+    else
+        return false;
+#else
+    // Else we just return true
+    return true;
+#endif
+}
+
+__inline__ __device__ bool cu_check_atom_type( int atom_type, int reference_type )
+{
+#ifdef SPIRIT_ENABLE_DEFECTS
+    // If defects are enabled we do a check if
+    //		atom types match.
+    if( atom_type == reference_type )
+        return true;
+    else
+        return false;
+#else
+    // Else we just return true
+    return true;
+#endif
+}
+
+// Calculates, for a spin i, a pair spin's index j.
+// This function takes into account boundary conditions and atom types and returns `-1` if any condition is not met.
+__inline__ __device__ int cu_idx_from_pair(
+    int ispin, const int * boundary_conditions, const int * n_cells, int N, const int * atom_types, const Pair & pair,
+    bool invert = false )
+{
+    // Invalid index if atom type of spin i is not correct
+    if( pair.i != ispin % N || !cu_check_atom_type( atom_types[ispin] ) )
+        return -1;
+
+    // Number of cells
+    auto & Na = n_cells[0];
+    auto & Nb = n_cells[1];
+    auto & Nc = n_cells[2];
+
+    // Invalid index if translations reach out over the lattice bounds
+    if( std::abs( pair.translations[0] ) > Na || std::abs( pair.translations[1] ) > Nb
+        || std::abs( pair.translations[2] ) > Nc )
+        return -1;
+
+    // Translations (cell) of spin i
+    int nic = ispin / ( N * Na * Nb );
+    int nib = ( ispin - nic * N * Na * Nb ) / ( N * Na );
+    int nia = ( ispin - nic * N * Na * Nb - nib * N * Na ) / N;
+
+    // Translations (cell) of spin j (possibly outside of non-periodical domain)
+    int pm = 1;
+    if( invert )
+        pm = -1;
+    int nja = nia + pm * pair.translations[0];
+    int njb = nib + pm * pair.translations[1];
+    int njc = nic + pm * pair.translations[2];
+
+    // Check boundary conditions: a
+    if( boundary_conditions[0] || ( 0 <= nja && nja < Na ) )
+    {
+        // Boundary conditions fulfilled
+        // Find the translations of spin j within the non-periodical domain
+        if( nja < 0 )
+            nja += Na;
+        // Calculate the correct index
+        if( nja >= Na )
+            nja -= Na;
+    }
+    else
+    {
+        // Boundary conditions not fulfilled
+        return -1;
+    }
+
+    // Check boundary conditions: b
+    if( boundary_conditions[1] || ( 0 <= njb && njb < Nb ) )
+    {
+        // Boundary conditions fulfilled
+        // Find the translations of spin j within the non-periodical domain
+        if( njb < 0 )
+            njb += Nb;
+        // Calculate the correct index
+        if( njb >= Nb )
+            njb -= Nb;
+    }
+    else
+    {
+        // Boundary conditions not fulfilled
+        return -1;
+    }
+
+    // Check boundary conditions: c
+    if( boundary_conditions[2] || ( 0 <= njc && njc < Nc ) )
+    {
+        // Boundary conditions fulfilled
+        // Find the translations of spin j within the non-periodical domain
+        if( njc < 0 )
+            njc += Nc;
+        // Calculate the correct index
+        if( njc >= Nc )
+            njc -= Nc;
+    }
+    else
+    {
+        // Boundary conditions not fulfilled
+        return -1;
+    }
+
+    // Calculate the index of spin j according to it's translations
+    int jspin = pair.j + (nja)*N + (njb)*N * Na + (njc)*N * Na * Nb;
+
+    if( !cu_check_atom_type( atom_types[jspin] ) )
+        return -1;
+
+    // Return a valid index
+    return jspin;
+}
+
+#endif
+
+inline std::array<int, 3> translations_from_idx( const intfield & n_cells, const int n_cell_atoms, int idx )
+{
+    std::array<int, 3> ret;
+    int Na = n_cells[0];
+    int Nb = n_cells[1];
+    int Nc = n_cells[2];
+    int N  = n_cell_atoms;
+
+    ret[2] = idx / ( N * Na * Nb );
+    ret[1] = ( idx - ret[2] * N * Na * Nb ) / ( N * Na );
+    ret[0] = ( idx - ret[2] * N * Na * Nb - ret[1] * N * Na ) / N;
+    return ret;
+}
+
+// Check atom types
+inline bool check_atom_type( int atom_type )
+{
+#ifdef SPIRIT_ENABLE_DEFECTS
+    // If defects are enabled we check for
+    //		vacancies (type < 0)
+    if( atom_type >= 0 )
+        return true;
+    else
+        return false;
+#else
+    // Else we just return true
+    return true;
+#endif
+}
+inline bool check_atom_type( int atom_type, int reference_type )
+{
+#ifdef SPIRIT_ENABLE_DEFECTS
+    // If defects are enabled we do a check if
+    //		atom types match.
+    if( atom_type == reference_type )
+        return true;
+    else
+        return false;
+#else
+    // Else we just return true
+    return true;
+#endif
+}
+
+// Calculates, for a spin i, a pair spin's index j.
+// This function takes into account boundary conditions and atom types and returns `-1` if any condition is not met.
+inline int idx_from_pair(
+    int ispin, const intfield & boundary_conditions, const intfield & n_cells, int N, const intfield & atom_types,
+    const Pair & pair, bool invert = false )
+{
+    // Invalid index if atom type of spin i is not correct
+    if( pair.i != ispin % N || !check_atom_type( atom_types[ispin] ) )
+        return -1;
+
+    // Number of cells
+    auto & Na = n_cells[0];
+    auto & Nb = n_cells[1];
+    auto & Nc = n_cells[2];
+
+    // Invalid index if translations reach out over the lattice bounds
+    if( std::abs( pair.translations[0] ) > Na || std::abs( pair.translations[1] ) > Nb
+        || std::abs( pair.translations[2] ) > Nc )
+        return -1;
+
+    // Translations (cell) of spin i
+    int nic = ispin / ( N * Na * Nb );
+    int nib = ( ispin - nic * N * Na * Nb ) / ( N * Na );
+    int nia = ( ispin - nic * N * Na * Nb - nib * N * Na ) / N;
+
+    int pm = 1;
+    if( invert )
+        pm = -1;
+    // Translations (cell) of spin j (possibly outside of non-periodical domain)
+    int nja = nia + pm * pair.translations[0];
+    int njb = nib + pm * pair.translations[1];
+    int njc = nic + pm * pair.translations[2];
+
+    // Check boundary conditions: a
+    if( boundary_conditions[0] || ( 0 <= nja && nja < Na ) )
+    {
+        // Boundary conditions fulfilled
+        // Find the translations of spin j within the non-periodical domain
+        if( nja < 0 )
+            nja += Na;
+        // Calculate the correct index
+        if( nja >= Na )
+            nja -= Na;
+    }
+    else
+    {
+        // Boundary conditions not fulfilled
+        return -1;
+    }
+
+    // Check boundary conditions: b
+    if( boundary_conditions[1] || ( 0 <= njb && njb < Nb ) )
+    {
+        // Boundary conditions fulfilled
+        // Find the translations of spin j within the non-periodical domain
+        if( njb < 0 )
+            njb += Nb;
+        // Calculate the correct index
+        if( njb >= Nb )
+            njb -= Nb;
+    }
+    else
+    {
+        // Boundary conditions not fulfilled
+        return -1;
+    }
+
+    // Check boundary conditions: c
+    if( boundary_conditions[2] || ( 0 <= njc && njc < Nc ) )
+    {
+        // Boundary conditions fulfilled
+        // Find the translations of spin j within the non-periodical domain
+        if( njc < 0 )
+            njc += Nc;
+        // Calculate the correct index
+        if( njc >= Nc )
+            njc -= Nc;
+    }
+    else
+    {
+        // Boundary conditions not fulfilled
+        return -1;
+    }
+
+    // Calculate the index of spin j according to it's translations
+    int jspin = pair.j + (nja)*N + (njb)*N * Na + (njc)*N * Na * Nb;
+
+    // Invalid index if atom type of spin j is not correct
+    if( !check_atom_type( atom_types[jspin] ) )
+        return -1;
+
+    // Return a valid index
+    return jspin;
+}
+
+}
+}
+
+#endif
diff --git a/core/include/engine/Vectormath.hpp b/core/include/engine/Vectormath.hpp
index 4b7c2dc80..38a541c6e 100644
--- a/core/include/engine/Vectormath.hpp
+++ b/core/include/engine/Vectormath.hpp
@@ -10,11 +10,13 @@
 #include <data/Geometry.hpp>
 #include <data/Spin_System.hpp>
 #include <engine/Vectormath_Defines.hpp>
+#include <engine/Indexing.hpp>
 
 namespace Engine
 {
 namespace Vectormath
 {
+
 // A "rotated view" into a vectorfield, with optional shifts applied before and after rotation.
 template<class vectorfield_t, class vector_t>
 class Rotated_View
@@ -54,479 +56,6 @@ void rotate( const vectorfield & v, const vectorfield & axis, const scalarfield
 // Decompose a vector into numbers of translations in a basis
 Vector3 decompose( const Vector3 & v, const std::vector<Vector3> & basis );
 
-/////////////////////////////////////////////////////////////////
-//////// Translating across the lattice
-
-// Note: translations must lie within bounds of n_cells
-inline int
-idx_from_translations( const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations )
-{
-    auto & Na = n_cells[0];
-    auto & Nb = n_cells[1];
-    auto & Nc = n_cells[2];
-    auto & N  = n_cell_atoms;
-
-    auto & da = translations[0];
-    auto & db = translations[1];
-    auto & dc = translations[2];
-
-    return da * N + db * N * Na + dc * N * Na * Nb;
-}
-
-#ifndef SPIRIT_USE_CUDA
-
-// Get the linear index in a n-D array where tupel contains the components in n-dimensions from fatest to slowest
-// varying and maxVal is the extent in every dimension
-inline int idx_from_tupel( const field<int> & tupel, const field<int> & maxVal )
-{
-    int idx  = 0;
-    int mult = 1;
-    for( int i = 0; i < tupel.size(); i++ )
-    {
-        idx += mult * tupel[i];
-        mult *= maxVal[i];
-    }
-    return idx;
-}
-
-// reverse of idx_from_tupel
-inline void tupel_from_idx( int & idx, field<int> & tupel, field<int> & maxVal )
-{
-    int idx_diff = idx;
-    int div      = 1;
-    for( int i = 0; i < tupel.size() - 1; i++ )
-        div *= maxVal[i];
-    for( int i = tupel.size() - 1; i > 0; i-- )
-    {
-        tupel[i] = idx_diff / div;
-        idx_diff -= tupel[i] * div;
-        div /= maxVal[i - 1];
-    }
-    tupel[0] = idx_diff / div;
-}
-
-inline int idx_from_translations(
-    const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations_i,
-    const std::array<int, 3> & translations )
-{
-    auto & Na = n_cells[0];
-    auto & Nb = n_cells[1];
-    auto & Nc = n_cells[2];
-    auto & N  = n_cell_atoms;
-
-    int da = translations_i[0] + translations[0];
-    int db = translations_i[1] + translations[1];
-    int dc = translations_i[2] + translations[2];
-
-    if( translations[0] < 0 )
-        da += N * Na;
-    if( translations[1] < 0 )
-        db += N * Na * Nb;
-    if( translations[2] < 0 )
-        dc += N * Na * Nb * Nc;
-
-    return ( da % Na ) * N + ( db % Nb ) * N * Na + ( dc % Nc ) * N * Na * Nb;
-}
-
-inline bool boundary_conditions_fulfilled(
-    const intfield & n_cells, const intfield & boundary_conditions, const std::array<int, 3> & translations_i,
-    const std::array<int, 3> & translations_j )
-{
-    int da = translations_i[0] + translations_j[0];
-    int db = translations_i[1] + translations_j[1];
-    int dc = translations_i[2] + translations_j[2];
-    return (
-        ( boundary_conditions[0] || ( 0 <= da && da < n_cells[0] ) )
-        && ( boundary_conditions[1] || ( 0 <= db && db < n_cells[1] ) )
-        && ( boundary_conditions[2] || ( 0 <= dc && dc < n_cells[2] ) ) );
-}
-
-// result = sum_i  f( vf1[i], vf2[i] )
-template<typename A, typename B, typename F>
-scalar reduce( const field<A> & vf1, const field<B> & vf2, const F & f )
-{
-    scalar res = 0;
-#pragma omp parallel for reduction( + : res )
-    for( unsigned int idx = 0; idx < vf1.size(); ++idx )
-    {
-        res += f( vf1[idx], vf2[idx] );
-    }
-    return res;
-}
-
-// vf1[i] = f( vf2[i] )
-template<typename A, typename B, typename F>
-void set( field<A> & vf1, const field<B> & vf2, const F & f )
-{
-#pragma omp parallel for
-    for( unsigned int idx = 0; idx < vf1.size(); ++idx )
-    {
-        vf1[idx] = f( vf2[idx] );
-    }
-}
-
-// f( vf1[idx], idx ) for all i
-template<typename F>
-void apply( int N, const F & f )
-{
-#pragma omp parallel for
-    for( unsigned int idx = 0; idx < N; ++idx )
-    {
-        f( idx );
-    }
-}
-
-#endif
-#ifdef SPIRIT_USE_CUDA
-
-// Get the linear index in a n-D array where tupel contains the components in n-dimensions from fatest to slowest
-// varying and maxVal is the extent in every dimension
-inline __device__ int cu_idx_from_tupel( field<int> & tupel, field<int> & maxVal )
-{
-    int idx  = 0;
-    int mult = 1;
-    for( int i = 0; i < tupel.size(); i++ )
-    {
-        idx += mult * tupel[i];
-        mult *= maxVal[i];
-    }
-    return idx;
-}
-
-// reverse of idx_from_tupel
-inline __device__ void cu_tupel_from_idx( int & idx, int * tupel, int * maxVal, int n )
-{
-    int idx_diff = idx;
-    int div      = 1;
-    for( int i = 0; i < n - 1; i++ )
-        div *= maxVal[i];
-    for( int i = n - 1; i > 0; i-- )
-    {
-        tupel[i] = idx_diff / div;
-        idx_diff -= tupel[i] * div;
-        div /= maxVal[i - 1];
-    }
-    tupel[0] = idx_diff / div;
-}
-
-inline void tupel_from_idx( int & idx, field<int> & tupel, const field<int> & maxVal )
-{
-    int idx_diff = idx;
-    int div      = 1;
-    for( int i = 0; i < maxVal.size() - 1; i++ )
-        div *= maxVal[i];
-    for( int i = maxVal.size() - 1; i > 0; i-- )
-    {
-        tupel[i] = idx_diff / div;
-        idx_diff -= tupel[i] * div;
-        div /= maxVal[i - 1];
-    }
-    tupel[0] = idx_diff / div;
-}
-
-inline int idx_from_translations(
-    const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations_i,
-    const int translations[3] )
-{
-    int Na = n_cells[0];
-    int Nb = n_cells[1];
-    int Nc = n_cells[2];
-    int N  = n_cell_atoms;
-
-    int da = translations_i[0] + translations[0];
-    int db = translations_i[1] + translations[1];
-    int dc = translations_i[2] + translations[2];
-
-    if( translations[0] < 0 )
-        da += N * Na;
-    if( translations[1] < 0 )
-        db += N * Na * Nb;
-    if( translations[2] < 0 )
-        dc += N * Na * Nb * Nc;
-
-    int idx = ( da % Na ) * N + ( db % Nb ) * N * Na + ( dc % Nc ) * N * Na * Nb;
-
-    return idx;
-}
-
-inline bool boundary_conditions_fulfilled(
-    const intfield & n_cells, const intfield & boundary_conditions, const std::array<int, 3> & translations_i,
-    const int translations_j[3] )
-{
-    int da = translations_i[0] + translations_j[0];
-    int db = translations_i[1] + translations_j[1];
-    int dc = translations_i[2] + translations_j[2];
-    return (
-        ( boundary_conditions[0] || ( 0 <= da && da < n_cells[0] ) )
-        && ( boundary_conditions[1] || ( 0 <= db && db < n_cells[1] ) )
-        && ( boundary_conditions[2] || ( 0 <= dc && dc < n_cells[2] ) ) );
-}
-
-__inline__ __device__ bool cu_check_atom_type( int atom_type )
-{
-#ifdef SPIRIT_ENABLE_DEFECTS
-    // If defects are enabled we check for
-    //		vacancies (type < 0)
-    if( atom_type >= 0 )
-        return true;
-    else
-        return false;
-#else
-    // Else we just return true
-    return true;
-#endif
-}
-
-__inline__ __device__ bool cu_check_atom_type( int atom_type, int reference_type )
-{
-#ifdef SPIRIT_ENABLE_DEFECTS
-    // If defects are enabled we do a check if
-    //		atom types match.
-    if( atom_type == reference_type )
-        return true;
-    else
-        return false;
-#else
-    // Else we just return true
-    return true;
-#endif
-}
-
-// Calculates, for a spin i, a pair spin's index j.
-// This function takes into account boundary conditions and atom types and returns `-1` if any condition is not met.
-__inline__ __device__ int cu_idx_from_pair(
-    int ispin, const int * boundary_conditions, const int * n_cells, int N, const int * atom_types, const Pair & pair,
-    bool invert = false )
-{
-    // Invalid index if atom type of spin i is not correct
-    if( pair.i != ispin % N || !cu_check_atom_type( atom_types[ispin] ) )
-        return -1;
-
-    // Number of cells
-    auto & Na = n_cells[0];
-    auto & Nb = n_cells[1];
-    auto & Nc = n_cells[2];
-
-    // Invalid index if translations reach out over the lattice bounds
-    if( std::abs( pair.translations[0] ) > Na || std::abs( pair.translations[1] ) > Nb
-        || std::abs( pair.translations[2] ) > Nc )
-        return -1;
-
-    // Translations (cell) of spin i
-    int nic = ispin / ( N * Na * Nb );
-    int nib = ( ispin - nic * N * Na * Nb ) / ( N * Na );
-    int nia = ( ispin - nic * N * Na * Nb - nib * N * Na ) / N;
-
-    // Translations (cell) of spin j (possibly outside of non-periodical domain)
-    int pm = 1;
-    if( invert )
-        pm = -1;
-    int nja = nia + pm * pair.translations[0];
-    int njb = nib + pm * pair.translations[1];
-    int njc = nic + pm * pair.translations[2];
-
-    // Check boundary conditions: a
-    if( boundary_conditions[0] || ( 0 <= nja && nja < Na ) )
-    {
-        // Boundary conditions fulfilled
-        // Find the translations of spin j within the non-periodical domain
-        if( nja < 0 )
-            nja += Na;
-        // Calculate the correct index
-        if( nja >= Na )
-            nja -= Na;
-    }
-    else
-    {
-        // Boundary conditions not fulfilled
-        return -1;
-    }
-
-    // Check boundary conditions: b
-    if( boundary_conditions[1] || ( 0 <= njb && njb < Nb ) )
-    {
-        // Boundary conditions fulfilled
-        // Find the translations of spin j within the non-periodical domain
-        if( njb < 0 )
-            njb += Nb;
-        // Calculate the correct index
-        if( njb >= Nb )
-            njb -= Nb;
-    }
-    else
-    {
-        // Boundary conditions not fulfilled
-        return -1;
-    }
-
-    // Check boundary conditions: c
-    if( boundary_conditions[2] || ( 0 <= njc && njc < Nc ) )
-    {
-        // Boundary conditions fulfilled
-        // Find the translations of spin j within the non-periodical domain
-        if( njc < 0 )
-            njc += Nc;
-        // Calculate the correct index
-        if( njc >= Nc )
-            njc -= Nc;
-    }
-    else
-    {
-        // Boundary conditions not fulfilled
-        return -1;
-    }
-
-    // Calculate the index of spin j according to it's translations
-    int jspin = pair.j + (nja)*N + (njb)*N * Na + (njc)*N * Na * Nb;
-
-    if( !cu_check_atom_type( atom_types[jspin] ) )
-        return -1;
-
-    // Return a valid index
-    return jspin;
-}
-
-#endif
-
-inline std::array<int, 3> translations_from_idx( const intfield & n_cells, const int n_cell_atoms, int idx )
-{
-    std::array<int, 3> ret;
-    int Na = n_cells[0];
-    int Nb = n_cells[1];
-    int Nc = n_cells[2];
-    int N  = n_cell_atoms;
-
-    ret[2] = idx / ( N * Na * Nb );
-    ret[1] = ( idx - ret[2] * N * Na * Nb ) / ( N * Na );
-    ret[0] = ( idx - ret[2] * N * Na * Nb - ret[1] * N * Na ) / N;
-    return ret;
-}
-
-// Check atom types
-inline bool check_atom_type( int atom_type )
-{
-#ifdef SPIRIT_ENABLE_DEFECTS
-    // If defects are enabled we check for
-    //		vacancies (type < 0)
-    if( atom_type >= 0 )
-        return true;
-    else
-        return false;
-#else
-    // Else we just return true
-    return true;
-#endif
-}
-inline bool check_atom_type( int atom_type, int reference_type )
-{
-#ifdef SPIRIT_ENABLE_DEFECTS
-    // If defects are enabled we do a check if
-    //		atom types match.
-    if( atom_type == reference_type )
-        return true;
-    else
-        return false;
-#else
-    // Else we just return true
-    return true;
-#endif
-}
-
-// Calculates, for a spin i, a pair spin's index j.
-// This function takes into account boundary conditions and atom types and returns `-1` if any condition is not met.
-inline int idx_from_pair(
-    int ispin, const intfield & boundary_conditions, const intfield & n_cells, int N, const intfield & atom_types,
-    const Pair & pair, bool invert = false )
-{
-    // Invalid index if atom type of spin i is not correct
-    if( pair.i != ispin % N || !check_atom_type( atom_types[ispin] ) )
-        return -1;
-
-    // Number of cells
-    auto & Na = n_cells[0];
-    auto & Nb = n_cells[1];
-    auto & Nc = n_cells[2];
-
-    // Invalid index if translations reach out over the lattice bounds
-    if( std::abs( pair.translations[0] ) > Na || std::abs( pair.translations[1] ) > Nb
-        || std::abs( pair.translations[2] ) > Nc )
-        return -1;
-
-    // Translations (cell) of spin i
-    int nic = ispin / ( N * Na * Nb );
-    int nib = ( ispin - nic * N * Na * Nb ) / ( N * Na );
-    int nia = ( ispin - nic * N * Na * Nb - nib * N * Na ) / N;
-
-    int pm = 1;
-    if( invert )
-        pm = -1;
-    // Translations (cell) of spin j (possibly outside of non-periodical domain)
-    int nja = nia + pm * pair.translations[0];
-    int njb = nib + pm * pair.translations[1];
-    int njc = nic + pm * pair.translations[2];
-
-    // Check boundary conditions: a
-    if( boundary_conditions[0] || ( 0 <= nja && nja < Na ) )
-    {
-        // Boundary conditions fulfilled
-        // Find the translations of spin j within the non-periodical domain
-        if( nja < 0 )
-            nja += Na;
-        // Calculate the correct index
-        if( nja >= Na )
-            nja -= Na;
-    }
-    else
-    {
-        // Boundary conditions not fulfilled
-        return -1;
-    }
-
-    // Check boundary conditions: b
-    if( boundary_conditions[1] || ( 0 <= njb && njb < Nb ) )
-    {
-        // Boundary conditions fulfilled
-        // Find the translations of spin j within the non-periodical domain
-        if( njb < 0 )
-            njb += Nb;
-        // Calculate the correct index
-        if( njb >= Nb )
-            njb -= Nb;
-    }
-    else
-    {
-        // Boundary conditions not fulfilled
-        return -1;
-    }
-
-    // Check boundary conditions: c
-    if( boundary_conditions[2] || ( 0 <= njc && njc < Nc ) )
-    {
-        // Boundary conditions fulfilled
-        // Find the translations of spin j within the non-periodical domain
-        if( njc < 0 )
-            njc += Nc;
-        // Calculate the correct index
-        if( njc >= Nc )
-            njc -= Nc;
-    }
-    else
-    {
-        // Boundary conditions not fulfilled
-        return -1;
-    }
-
-    // Calculate the index of spin j according to it's translations
-    int jspin = pair.j + (nja)*N + (njb)*N * Na + (njc)*N * Na * Nb;
-
-    // Invalid index if atom type of spin j is not correct
-    if( !check_atom_type( atom_types[jspin] ) )
-        return -1;
-
-    // Return a valid index
-    return jspin;
-}
-
 /////////////////////////////////////////////////////////////////
 //////// Vectorfield Math - special stuff