diff --git a/distopia/__init__.py b/distopia/__init__.py
index 6e07abf3..8a062d6d 100644
--- a/distopia/__init__.py
+++ b/distopia/__init__.py
@@ -16,6 +16,9 @@
     calc_self_distance_array_no_box,
     calc_self_distance_array_ortho,
     calc_self_distance_array_triclinic,
+    calc_bonds_no_box_idx,
+    calc_bonds_ortho_idx,
+    calc_bonds_triclinic_idx,
 )
 
 __all__ = [
@@ -34,5 +37,8 @@
     'calc_self_distance_array_no_box',
     'calc_self_distance_array_ortho',
     'calc_self_distance_array_triclinic',
+    'calc_bonds_no_box_idx',
+    'calc_bonds_ortho_idx',
+    'calc_bonds_triclinic_idx',
     '__version__',
 ]
\ No newline at end of file
diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx
index 5b0d077b..78519701 100644
--- a/distopia/_distopia.pyx
+++ b/distopia/_distopia.pyx
@@ -13,7 +13,7 @@ The python functions for distopia
 import numpy as np
 cimport cython
 cimport numpy as cnp
-from cython cimport floating
+from cython cimport floating, integral
 cnp.import_array()
 
 cdef extern from "distopia.h" namespace "distopia" nogil:
@@ -129,6 +129,95 @@ cdef extern from "distopia.h" namespace "distopia" nogil:
         const T *box,
         T *out,
     )
+    void CalcBondsNoBoxIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        size_t n,
+        T *out,
+    )
+
+    void CalcBondsOrthoIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcBondsTriclinicIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcAnglesNoBoxIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        size_t n,
+        T *out,
+    )
+
+    void CalcAnglesOrthoIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcAnglesTriclinicIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcDihedralsNoBoxIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        const int* d_idx,
+        size_t n,
+        T *out,
+    )
+
+    void CalcDihedralsOrthoIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        const int* d_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+
+    void CalcDihedralsTriclinicIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        const int* d_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+
 
 def get_n_float_lanes():
     """The number of floats per register distopia will handle on this system"""
@@ -890,4 +979,152 @@ def calc_self_distance_array_triclinic(
                                &box[0][0],
                                &results_view[0])
 
-    return np.array(results).reshape(coords0.shape[0], coords0.shape[0])
\ No newline at end of file
+    return np.array(results).reshape(coords0.shape[0], coords0.shape[0])
+
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def calc_bonds_no_box_idx(
+     floating[:, ::1] coords,
+     int[::1] a_idx,
+     int[::1] b_idx,
+     floating[::1] results=None):
+    """Calculate pairwise distances between coords[a_idx] and coords[b_idx] with no periodic boundary conditions
+
+    Parameters
+    ----------
+    coords : float32 or float64 array
+      must be same length and dtype
+    a_idx, b_idx : int array
+      must be same length and dtype
+    results: float32 or float64 array (optional)
+      array to store results in, must be same size and dtype as a_idx/b_idx
+
+    Returns
+    -------
+    distances : float32 or float64 array
+      same length and dtype as a_idx/b_idx
+    """
+    cdef floating[::1] results_view
+    cdef size_t nvals = a_idx.shape[0]
+    cdef cnp.npy_intp[1] dims
+    dims[0] = <ssize_t > nvals  # FIXME truncation?
+
+    _check_shapes(a_idx, b_idx)
+
+    if results is None:
+        if floating is float:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
+        else:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0)
+
+    else:
+        _check_results(results, nvals)
+
+    results_view = results
+
+    CalcBondsNoBoxIdx(& coords[0][0], & a_idx[0], & b_idx[0], nvals, & results_view[0])
+
+    return np.array(results)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def calc_bonds_ortho_idx(
+     floating[:, ::1] coords,
+     int[::1] a_idx,
+     int[::1] b_idx,
+     floating[::1] box,
+     floating[::1] results=None):
+    """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under orthorhombic boundary conditions
+
+    Parameters
+    ----------
+    coords : float32 or float64 array
+      must be same length and dtype
+    a_idx, b_idx : int array
+      must be same length and dtype
+    box : float32 or float64 array
+      orthorhombic periodic boundary dimensions in [L, L, L] format
+    results: float32 or float64 array (optional)
+      array to store results in, must be same size and dtype as a_idx/b_idx
+
+    Returns
+    -------
+    distances : float32 or float64 array
+      same length and dtype as a_idx/b_idx
+    """
+    cdef floating[::1] results_view
+    cdef size_t nvals = a_idx.shape[0]
+    cdef cnp.npy_intp[1] dims
+    dims[0] = <ssize_t > nvals  # FIXME truncation?
+
+    _check_shapes(a_idx, b_idx)
+
+
+    if results is None:
+        if floating is float:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
+        else:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0)
+
+    else:
+        _check_results(results, nvals)
+
+    results_view = results
+
+    CalcBondsOrthoIdx(& coords[0][0], & a_idx[0], & b_idx[0], nvals, & box[0], & results_view[0])
+
+    return np.array(results)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def calc_bonds_triclinic_idx(
+     floating[:, ::1] coords,
+     int[::1] a_idx,
+     int[::1] b_idx,
+     floating[:, ::1] box,
+     floating[::1] results=None):
+    """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under triclinic boundary conditions
+
+    Parameters
+    ----------
+    coords : float32 or float64 array
+      must be same length and dtype
+    a_idx, b_idx : int array
+      must be same length and dtype
+    box : float32 or float64 array
+      periodic boundary dimensions, in 3x3 format
+    results: float32 or float64 array (optional)
+      array to store results in, must be same size and dtype as a_idx/b_idx
+
+    Returns
+    -------
+    distances : float32 or float64 array
+      same length and dtype as a_idx/b_idx
+    """
+    cdef floating[::1] results_view
+    cdef size_t nvals = a_idx.shape[0]
+    cdef cnp.npy_intp[1] dims
+    dims[0] = <ssize_t > nvals  # FIXME truncation?
+
+    _check_shapes(a_idx, b_idx)
+
+
+    if results is None:
+        if floating is float:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
+        else:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0)
+
+    else:
+        _check_results(results, nvals)
+
+    results_view = results
+
+    CalcBondsTriclinicIdx(& coords[0][0], & a_idx[0], & b_idx[0], nvals, & box[0][0], & results_view[0])
+
+    return np.array(results)
+
diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py
index 798780d8..11683d24 100644
--- a/distopia/tests/test_distopia.py
+++ b/distopia/tests/test_distopia.py
@@ -240,6 +240,36 @@ def test_triclinic_bad_result(self):
             distopia.calc_self_distance_array_triclinic(c0, box, results=np.empty(1, dtype=np.float32))
 
 
+
+class TestDistancesIdx:
+    def arange_input(self, N, dtype):
+        return np.arange(3 * N, dtype=dtype).reshape(N, 3)
+
+    def result_shim(self, make_arr, N, dtype):
+        if not make_arr:
+            return None
+        else:
+            return np.empty(N, dtype=dtype)
+
+
+    @pytest.mark.parametrize("dtype", (np.float32, np.float64))
+    @pytest.mark.parametrize("N", (3, 6, 10, 1000, 10000))
+    @pytest.mark.parametrize("use_result_buffer", (True, False))
+    def test_calc_bonds_no_box_idx(self, N, use_result_buffer, dtype):
+        c = self.arange_input(N, dtype)
+
+        a_idx = np.asarray(np.arange(N -3)).astype(np.int32)
+        b_idx = np.asarray(np.arange(N -3)).astype(np.int32)
+        result_buffer = self.result_shim(use_result_buffer, len(a_idx), dtype)
+        result = distopia.calc_bonds_no_box_idx(
+            c, a_idx, b_idx, results=result_buffer
+        )
+        assert_allclose(result, np.zeros(N))
+        # check dtype of result
+        assert result.dtype == dtype
+    
+
+
 class TestMDA:
     """
     Copy of some of the tests from MDAnalysisTests repo
@@ -381,4 +411,8 @@ def test_periodic_angles(self, box_bonds, positions_angles, dtype):
         test3 = distopia.calc_angles_ortho(a, b, c2, box=box)
 
         for val in [test1, test2, test3]:
-            assert_almost_equal(ref, val, self.prec, err_msg="Min image in angle calculation failed")
\ No newline at end of file
+            assert_almost_equal(ref, val, self.prec, err_msg="Min image in angle calculation failed")
+
+
+
+
diff --git a/libdistopia/CMakeLists.txt b/libdistopia/CMakeLists.txt
index 23b9a45b..ef05c6d0 100644
--- a/libdistopia/CMakeLists.txt
+++ b/libdistopia/CMakeLists.txt
@@ -41,10 +41,14 @@ Include(GoogleTest)
 add_subdirectory("googletest")
 enable_testing()
 include_directories(${gtest_SOURCE_DIR}/include ${gtest_SOURCE_DIR})
+include_directories(${gmock_SOURCE_DIR}/include ${gmock_SOURCE_DIR})
+
 
 add_executable(test)
 target_sources(test PRIVATE "test/test.cpp")
 target_link_libraries(test PUBLIC gtest gtest_main)
+target_link_libraries(test PUBLIC gmock gmock_main)
+
 target_link_libraries(test PUBLIC libdistopia)
 target_include_directories(test PUBLIC  ${CMAKE_SOURCE_DIR})
 gtest_discover_tests(test WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY})
@@ -58,6 +62,7 @@ target_include_directories(targets PUBLIC  ${CMAKE_SOURCE_DIR})
 add_executable(test_mda_match)
 target_sources(test_mda_match PRIVATE "test/test_mda_match.cpp")
 target_link_libraries(test_mda_match PUBLIC gtest gtest_main)
+target_link_libraries(test PUBLIC gmock gmock_main)
 target_link_libraries(test_mda_match PUBLIC libdistopia)
 target_include_directories(test_mda_match PUBLIC  ${CMAKE_SOURCE_DIR})
 gtest_discover_tests(test_mda_match WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY})
diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp
index b1187fa0..212bfb22 100644
--- a/libdistopia/src/distopia.cpp
+++ b/libdistopia/src/distopia.cpp
@@ -54,7 +54,7 @@ namespace distopia {
                 acc = hn::MulAdd(dx, dx, acc);
                 acc = hn::MulAdd(dy, dy, acc);
                 acc = hn::MulAdd(dz, dz, acc);
-
+                hn::Print(d, "rsq", acc);
                 return hn::Sqrt(acc);
             }
         };
@@ -492,7 +492,7 @@ namespace distopia {
 
             y = y / vb_norm;
 
-            return   hn::Neg(hn::Atan2(d, y, x));
+            return hn::Neg(hn::Atan2(d, y, x));
         }
 
         template <typename T, typename B>
@@ -1607,74 +1607,132 @@ namespace distopia {
     }
     HWY_DLLEXPORT template <> void CalcBondsNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx,
                                                      int n, float *out) {
+
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcBondsNoBoxIdxSingle(coords, a_idx, b_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsNoBoxIdxSingle)(coords, a_idx, b_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx,
                                                      int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcBondsNoBoxIdxDouble(coords, a_idx, b_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsNoBoxIdxDouble)(coords, a_idx, b_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsOrthoIdx(const float *coords, const int *a_idx, const int *b_idx,
                                                      int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcBondsOrthoIdxSingle(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsOrthoIdxSingle)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsOrthoIdx(const double *coords, const int *a_idx, const int *b_idx,
                                                      int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcBondsOrthoIdxDouble(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsOrthoIdxDouble)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsTriclinicIdx(const float *coords, const int *a_idx, const int *b_idx,
                                                          int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcBondsTriclinicIdxSingle(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsTriclinicIdxSingle)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx,
                                                          int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcBondsTriclinicIdxDouble(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsTriclinicIdxDouble)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcAnglesNoBoxIdxSingle(coords, a_idx, b_idx, c_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesNoBoxIdxSingle)(coords, a_idx, b_idx, c_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesNoBoxIdxDouble(coords, a_idx, b_idx, c_idx, n, out);
+        }
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesNoBoxIdxDouble(coords, a_idx, b_idx, c_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesNoBoxIdxDouble)(coords, a_idx, b_idx, c_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesOrthoIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcAnglesOrthoIdxSingle(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesOrthoIdxSingle)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesOrthoIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesOrthoIdxDouble(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesOrthoIdxDouble)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesTriclinicIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                           int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcAnglesTriclinicIdxSingle(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesTriclinicIdxSingle)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                           int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesTriclinicIdxDouble(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesTriclinicIdxDouble)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsNoBoxIdxSingle(coords, a_idx, b_idx, c_idx, d_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsNoBoxIdxSingle)(coords, a_idx, b_idx, c_idx, d_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsNoBoxIdxDouble(coords, a_idx, b_idx, c_idx, d_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsNoBoxIdxDouble)(coords, a_idx, b_idx, c_idx, d_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsOrthoIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsOrthoIdxSingle(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsOrthoIdxSingle)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsOrthoIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsOrthoIdxDouble(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsOrthoIdxDouble)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsTriclinicIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                              int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsTriclinicIdxSingle(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsTriclinicIdxSingle)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                              int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsTriclinicIdxDouble(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsTriclinicIdxDouble)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx, int na, int nb, float *out) {
@@ -1684,7 +1742,7 @@ namespace distopia {
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayNoBoxIdxSingle)(coords, a_idx, b_idx, na, nb, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx, int na, int nb, double *out) {
-        if (nb < GetNFloatLanes()) {
+        if (nb < GetNDoubleLanes()) {
             return distopia::N_SCALAR::CalcDistanceArrayNoBoxIdxDouble(coords, a_idx, b_idx, na, nb, out);
         }
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayNoBoxIdxDouble)(coords, a_idx, b_idx, na, nb, out);
@@ -1696,7 +1754,7 @@ namespace distopia {
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayOrthoIdxSingle)(coords, a_idx, b_idx, na, nb, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayOrthoIdx(const double *coords, const int *a_idx, const int *b_idx, int na, int nb, const double *box, double *out) {
-        if (nb < GetNFloatLanes()) {
+        if (nb < GetNDoubleLanes()) {
             return distopia::N_SCALAR::CalcDistanceArrayOrthoIdxDouble(coords, a_idx, b_idx, na, nb, box, out);
         }
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayOrthoIdxDouble)(coords, a_idx, b_idx, na, nb, box, out);
@@ -1708,27 +1766,45 @@ namespace distopia {
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayTriclinicIdxSingle)(coords, a_idx, b_idx, na, nb, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx, int na, int nb, const double *box, double *out) {
-        if (nb < GetNFloatLanes()) {
+        if (nb < GetNDoubleLanes()) {
             return distopia::N_SCALAR::CalcDistanceArrayTriclinicIdxDouble(coords, a_idx, b_idx, na, nb, box, out);
         }
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayTriclinicIdxDouble)(coords, a_idx, b_idx, na, nb, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayNoBoxIdx(const float *coords, const int *idx, int n, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayNoBoxIdxSingle(coords, idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayNoBoxIdxSingle)(coords, idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayNoBoxIdx(const double *coords, const int *idx, int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayNoBoxIdxDouble(coords, idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayNoBoxIdxDouble)(coords, idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayOrthoIdx(const float *coords, const int *idx, int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayOrthoIdxSingle(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayOrthoIdxSingle)(coords, idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayOrthoIdx(const double *coords, const int *idx, int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayOrthoIdxDouble(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayOrthoIdxDouble)(coords, idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayTriclinicIdx(const float *coords, const int *idx, int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayTriclinicIdxSingle(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayTriclinicIdxSingle)(coords, idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayTriclinicIdx(const double *coords, const int *idx, int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayTriclinicIdxDouble(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayTriclinicIdxDouble)(coords, idx, n, box, out);
     }
 
diff --git a/libdistopia/test/compare/calc_distances.h b/libdistopia/test/compare/calc_distances.h
index 4be27d57..330b0e30 100644
--- a/libdistopia/test/compare/calc_distances.h
+++ b/libdistopia/test/compare/calc_distances.h
@@ -607,7 +607,10 @@ static void _calc_bond_distance(ScalarToCoordinateT<T>* atom1, ScalarToCoordinat
     dx[1] = atom1[i][1] - atom2[i][1];
     dx[2] = atom1[i][2] - atom2[i][2];
     T rsq = (dx[0]*dx[0])+(dx[1]*dx[1])+(dx[2]*dx[2]);
+    std::cout << "rsq = " << rsq << std::endl;
+
     *(distances+i) = sqrt(rsq);
+    std::cout << "distances[" << i << "] = " << sqrt(rsq) << std::endl;
   }
 }
 
diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp
index e226ec6c..c71cd234 100644
--- a/libdistopia/test/test_mda_match.cpp
+++ b/libdistopia/test/test_mda_match.cpp
@@ -1,6 +1,7 @@
 #include <cmath>
 #include <iostream>
 #include <numeric>
+#include <stdlib.h>     
 
 #include "gtest/gtest.h"
 #include "distopia.h"
@@ -14,8 +15,8 @@ typedef Types<float, double> ScalarTypes;
 
 // constants
 constexpr int BOXSIZE = 30;
-constexpr int NRESULTS = 10;
-constexpr int NINDICIES = 5;
+constexpr int NRESULTS = 250;
+constexpr int NINDICIES = 100;
 constexpr double abs_err = 1.0e-4;
 
 
@@ -39,7 +40,7 @@ TYPED_TEST(CoordinatesTest, CalcBondsOrthoMatchesMDA)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -56,7 +57,7 @@ TYPED_TEST(CoordinatesTest, CalcBondsNoBoxMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -92,7 +93,7 @@ TYPED_TEST(CoordinatesTest, CalcBondsTriclinicMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -111,7 +112,7 @@ TYPED_TEST(CoordinatesTest, CalcAnglesOrthoMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -128,7 +129,7 @@ TYPED_TEST(CoordinatesTest, CalcAnglesNoBoxMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -157,7 +158,7 @@ TYPED_TEST(CoordinatesTest, CalcAnglesTriclinicMatchesMDA)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -174,7 +175,7 @@ TYPED_TEST(CoordinatesTest, CalcDihedralsOrthoMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -191,7 +192,7 @@ TYPED_TEST(CoordinatesTest, CalcDihedralsNoBoxMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -221,7 +222,7 @@ TYPED_TEST(CoordinatesTest, CalcDihedralsTriclinicMatchesMDA)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -241,7 +242,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDA) {
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 
 }
@@ -261,7 +262,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDAWide) {
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 
 }
@@ -282,7 +283,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDAScalarPath)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -300,7 +301,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayNoBoxMatchesMDA) {
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -320,7 +321,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayNoBoxMatchesMDAScalarPath)
 
     for (std::size_t i = 0; i <  this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -338,7 +339,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayTriclinicMatchesMDA) {
 
     for (std::size_t i = 0;  i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -358,7 +359,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayTriclinicMatchesMDAScalarP
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -377,7 +378,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayNoBox) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -397,7 +398,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayOrtho) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -417,7 +418,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayTriclinic) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -437,7 +438,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayNoBoxScalar) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -457,7 +458,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayOrthoScalar) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -477,7 +478,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayTriclinicScalar) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -492,13 +493,40 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcBondsNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
+    // distopia::CalcBondsNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
+    using ctype = ScalarToCoordinateT<TypeParam>;
+
+    _calc_bond_distance((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->ref_results);
 
     // test the idx
     distopia::CalcBondsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+    }
+
+    // print 
+    for (std::size_t i=0; i<this->nidx; i++) {
+        std::cout << "i: " << i << " res: " << this->results[i] << " ressq: " <<  std::pow(this->results[i], 2) << " ref: " << this->ref_results[i] << " refsq: " << std::pow(this->ref_results[i], 2) << " near: " << isnear(this->results[i], this->ref_results[i], abs_err) << std::endl;
+    }
+
+}
+
+TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
+    int ncoords = 250;
+    // 3 is always a wierd remainder case due to 
+    int nidx = 3;
+
+    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+    using ctype = ScalarToCoordinateT<TypeParam>;
+
+    // reference result
+    _calc_bond_distance((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->ref_results);
+    // test the idx
+    distopia::CalcBondsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
+
+    for (std::size_t i=0; i<this->nidx; i++) {
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -508,33 +536,35 @@ TYPED_TEST(CoordinatesIdx, CalcBondsOrthoIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+    using ctype = ScalarToCoordinateT<TypeParam>;
 
     // reference result
-    distopia::CalcBondsOrtho(this->a_coords_contig, this->b_coords_contig, this->nidx, this->box, this->ref_results);
+    _calc_bond_distance_ortho((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->box, this->ref_results);
 
     // test the idx
-    distopia::CalcBondsOrthoIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->box, this->results);
+    distopia::CalcBondsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
 
+
 TYPED_TEST(CoordinatesIdx, CalcBondsTriclinicIdx) {
     int ncoords = 250;
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+    using ctype = ScalarToCoordinateT<TypeParam>;
 
     // reference result
-    distopia::CalcBondsTriclinic(this->a_coords_contig, this->b_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
-
+    _calc_bond_distance_triclinic((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
     // test the idx
-    distopia::CalcBondsTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcBondsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -544,15 +574,14 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesNoBoxIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcAnglesNoBox(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->ref_results);
-
+    _calc_angle((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, this->nidx, this->ref_results);
     // test the idx
-    distopia::CalcAnglesNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->results);
+    distopia::CalcAnglesNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -562,15 +591,14 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesOrthoIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcAnglesOrtho(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->box, this->ref_results);
-
+    _calc_angle_ortho((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, this->nidx, this->box, this->ref_results);
     // test the idx
-    distopia::CalcAnglesOrthoIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->box, this->results);
+    distopia::CalcAnglesOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -578,17 +606,17 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesOrthoIdx) {
 TYPED_TEST(CoordinatesIdx, CalcAnglesTriclinicIdx) {
     int ncoords = 250;
     int nidx = 100;
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcAnglesTriclinic(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+    _calc_angle_triclinic((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
     // test the idx
-    distopia::CalcAnglesTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcAnglesTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -598,15 +626,14 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsNoBoxIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcDihedralsNoBox(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->ref_results);
-
+    _calc_dihedral((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, (ctype*)this->d_coords_contig, this->nidx, this->ref_results);
     // test the idx
-    distopia::CalcDihedralsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->results);
+    distopia::CalcDihedralsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -616,15 +643,14 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsOrthoIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcDihedralsOrtho(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->box, this->ref_results);
-
+    _calc_dihedral_ortho((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, (ctype*)this->d_coords_contig, this->nidx, this->box, this->ref_results);
     // test the idx
-    distopia::CalcDihedralsOrthoIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->box, this->results);
+    distopia::CalcDihedralsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -634,134 +660,133 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsTriclinicIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcDihedralsTriclinic(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
-
+    _calc_dihedral_triclinic((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, (ctype*)this->d_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
     // test the idx
-    distopia::CalcDihedralsTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcDihedralsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
 
-TYPED_TEST(CoordinatesIdx, DistanceArrayNoBoxIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, DistanceArrayNoBoxIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcDistanceArrayNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->ref_results);
+//     // ref
+//     distopia::CalcDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->ref_results);
 
-    // test
-    distopia::CalcDistanceArrayNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->results);
+//     // test
+//     distopia::CalcDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->results);
 
-    size_t n_results = nidx * nidx;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
+//     size_t n_results = nidx * nidx;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
 
-}
+// }
 
 
-TYPED_TEST(CoordinatesIdx, DistanceArrayOrthoIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, DistanceArrayOrthoIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcDistanceArrayOrtho(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->box, this->ref_results);
+//     // ref
+//     distopia::CalcDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->box, this->ref_results);
 
-    // test
-    distopia::CalcDistanceArrayOrthoIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->box, this->results);
+//     // test
+//     distopia::CalcDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->box, this->results);
 
-    size_t n_results = nidx * nidx;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
+//     size_t n_results = nidx * nidx;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
 
-}
+// }
 
 
-TYPED_TEST(CoordinatesIdx, DistanceArrayTriclinicIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, DistanceArrayTriclinicIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcDistanceArrayTriclinic(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->triclinic_box, this->ref_results);
+//     // ref
+//     distopia::CalcDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->triclinic_box, this->ref_results);
 
-    // test
-    distopia::CalcDistanceArrayTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->triclinic_box, this->results);
+//     // test
+//     distopia::CalcDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->triclinic_box, this->results);
 
-    size_t n_results = nidx * nidx;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
+//     size_t n_results = nidx * nidx;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
 
-}
+// }
 
 
-TYPED_TEST(CoordinatesIdx, SelfDistanceArrayNoBoxIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, SelfDistanceArrayNoBoxIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcSelfDistanceArrayNoBox(this->a_coords_contig, this->nidx, this->ref_results);
+//     // ref
+//     distopia::CalcSelfDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->nidx, this->ref_results);
 
-    // test
-    distopia::CalcSelfDistanceArrayNoBoxIdx(this->coords, this->a_idx, this->nidx, this->results);
+//     // test
+//     distopia::CalcSelfDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->results);
 
 
-    size_t n_results = nidx * (nidx-1) / 2;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
-}
+//     size_t n_results = nidx * (nidx-1) / 2;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
+// }
 
 
-TYPED_TEST(CoordinatesIdx, SelfDistanceArrayOrthoIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, SelfDistanceArrayOrthoIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcSelfDistanceArrayOrtho(this->a_coords_contig, this->nidx, this->box, this->ref_results);
+//     // ref
+//     distopia::CalcSelfDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->nidx, this->box, this->ref_results);
 
-    // test
-    distopia::CalcSelfDistanceArrayOrthoIdx(this->coords, this->a_idx, this->nidx, this->box, this->results);
+//     // test
+//     distopia::CalcSelfDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->box, this->results);
 
 
-    size_t n_results = nidx * (nidx-1) / 2;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
-}
+//     size_t n_results = nidx * (nidx-1) / 2;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
+// }
 
 
-TYPED_TEST(CoordinatesIdx, SelfDistanceArrayTriclinicIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, SelfDistanceArrayTriclinicIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcSelfDistanceArrayTriclinic(this->a_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+//     // ref
+//     distopia::CalcSelfDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
-    // test
-    distopia::CalcSelfDistanceArrayTriclinicIdx(this->coords, this->a_idx, this->nidx, this->triclinic_box, this->results);
+//     // test
+//     distopia::CalcSelfDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->triclinic_box, this->results);
 
 
-    size_t n_results = nidx * (nidx-1) / 2;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
-}
\ No newline at end of file
+//     size_t n_results = nidx * (nidx-1) / 2;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
+// }
\ No newline at end of file
diff --git a/libdistopia/test/test_utils.h b/libdistopia/test/test_utils.h
index 469da52f..2057686b 100644
--- a/libdistopia/test/test_utils.h
+++ b/libdistopia/test/test_utils.h
@@ -4,6 +4,9 @@
 #include <random>
 
 #include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+
 // creates nrandom floating points between pos and neg limit
 template <typename T>
 void RandomFloatingPoint(T *target, const int nrandom, const int neglimit,
@@ -33,21 +36,33 @@ void RandomInt(std::size_t *target, const int nrandom, const int neglimit,
 
 inline void EXPECT_SCALAR_EQ(float result, float ref)
 {
-    EXPECT_FLOAT_EQ(result, ref);
+
+    EXPECT_THAT(result, ::testing::NanSensitiveFloatEq(ref));
 }
 
 inline void EXPECT_SCALAR_EQ(double result, double ref)
 {
-    EXPECT_DOUBLE_EQ(result, ref);
+
+    EXPECT_THAT(result, ::testing::NanSensitiveDoubleEq(ref));
 }
 
 inline void EXPECT_SCALAR_NEAR(float result, float ref, float tol)
 {
-    EXPECT_NEAR(result, ref, tol);
+
+    EXPECT_THAT(result, ::testing::NanSensitiveFloatNear(ref, tol));
 }
 inline void EXPECT_SCALAR_NEAR(double result, double ref, float tol)
 {
-    EXPECT_NEAR(result, ref, tol);
+
+    EXPECT_THAT(result, ::testing::NanSensitiveDoubleNear(ref, tol));
+}
+
+
+// isnear 
+
+bool isnear(float a, float b, float tol)
+{
+    return std::abs(a - b) < tol;
 }
 
 
@@ -70,4 +85,6 @@ void pretty_print_matrix(T *matrix, int rows, int cols)
 
 
 
+
+
 #endif // DISTOPIA_TEST_UTILS_H
\ No newline at end of file