diff --git a/.github/workflows/cmake.yml b/.github/workflows/build_and_test.yml
similarity index 95%
rename from .github/workflows/cmake.yml
rename to .github/workflows/build_and_test.yml
index 6f463f2e..b89d1410 100644
--- a/.github/workflows/cmake.yml
+++ b/.github/workflows/build_and_test.yml
@@ -1,4 +1,4 @@
-name: CMake
+name: build_and_test
 
 on:
   push:
@@ -79,7 +79,7 @@ jobs:
       # Execute tests defined by the CMake configuration.  
       # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
       run:  ctest --test-dir _skbuild/*/cmake-build/libdistopia
-
+    
   pip-install:
     # A pure Python install, which relies purely on pyproject.toml contents
     runs-on: ${{ matrix.os }}
@@ -109,3 +109,10 @@ jobs:
 
     - name: test
       run: ctest --test-dir _skbuild/*/cmake-build/libdistopia
+
+    - name: install_pytest
+      run: pip install pytest
+
+    - name: python_test
+      # run python API tests
+      run: pytest distopia/tests
diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py
new file mode 100644
index 00000000..08bb39da
--- /dev/null
+++ b/distopia/tests/test_distopia.py
@@ -0,0 +1,94 @@
+import pytest
+import distopia
+import numpy as np
+from numpy.testing import assert_allclose
+
+
+"""
+Majority of detailed testing is done at the C++ level.
+This is primarily to make sure that the Python API works as expected.
+"""
+
+
+class TestDistances:
+
+    def arange_input(self, N, dtype):
+        return np.arange(3*N, dtype=dtype).reshape(N,3)
+
+    def idx_input(self, N):
+        idxs = np.zeros((N, 2), dtype=np.ulonglong)
+        for i in range(N):
+            idxs[i,0]=i
+            idxs[i,1]=i+N
+        return np.ravel(idxs)
+
+    @pytest.mark.parametrize('box', ([10, 10, 10], [100, 20, 10]))
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_ortho_float_all_zero(self, N, box):
+        c0 = self.arange_input(N, np.float32)
+        c1 = self.arange_input(N, np.float32)
+        result = distopia.calc_bonds_ortho_float(
+            c0, c1, np.asarray(box, dtype=np.float32))
+        assert_allclose(result, np.zeros(N))
+
+    @pytest.mark.parametrize('box', ([10, 10, 10], [100, 20, 10]))
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_ortho_double_all_zero(self, N, box):
+        c0 = self.arange_input(N, np.float64)
+        c1 = self.arange_input(N, np.float64)
+        result = distopia.calc_bonds_ortho_double(
+            c0, c1, np.asarray(box, dtype=np.float64))
+        assert_allclose(result, np.zeros(N))
+
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_nobox_float_all_zero(self, N):
+        c0 = self.arange_input(N, np.float32)
+        c1 = self.arange_input(N, np.float32)
+        result = distopia.calc_bonds_no_box_float(c0, c1)
+        assert_allclose(result, np.zeros(N))
+    
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_nobox_double_all_zero(self, N):
+        c0 = self.arange_input(N, np.float64)
+        c1 = self.arange_input(N, np.float64)
+        result = distopia.calc_bonds_no_box_double(c0, c1)
+        assert_allclose(result, np.zeros(N))
+
+    @pytest.mark.parametrize('box', ([10, 10, 10], [100, 20, 10]))
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_idx_ortho_float_all_zero(self, N, box):
+        c0 = self.arange_input(N, np.float32)
+        c1 = self.arange_input(N, np.float32)
+        idx = self.idx_input(N)
+        coords = np.vstack((c0, c1))
+        result = distopia.calc_bonds_idx_ortho_float(coords, idx, np.asarray(box).astype(np.float32))
+        assert_allclose(result, np.zeros(N))
+
+    @pytest.mark.parametrize('box', ([10, 10, 10], [100, 20, 10]))
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_idx_ortho_double_all_zero(self, N, box):
+        c0 = self.arange_input(N, np.float64)
+        c1 = self.arange_input(N, np.float64)
+        idx = self.idx_input(N)
+        coords = np.vstack((c0, c1))
+        result = distopia.calc_bonds_idx_ortho_double(coords, idx, np.asarray(box).astype(np.float64))
+        assert_allclose(result, np.zeros(N))
+    
+
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_idx_no_box_float_all_zero(self, N):
+        c0 = self.arange_input(N, np.float32)
+        c1 = self.arange_input(N, np.float32)
+        idx = self.idx_input(N)
+        coords = np.vstack((c0, c1))
+        result = distopia.calc_bonds_idx_no_box_float(coords, idx)
+        assert_allclose(result, np.zeros(N))
+
+    @pytest.mark.parametrize('N', (0, 10, 1000, 10000))
+    def test_calc_bonds_idx_no_box_double_all_zero(self, N):
+        c0 = self.arange_input(N, np.float64)
+        c1 = self.arange_input(N, np.float64)
+        idx = self.idx_input(N)
+        coords = np.vstack((c0, c1))
+        result = distopia.calc_bonds_idx_no_box_double(coords, idx)
+        assert_allclose(result, np.zeros(N))
diff --git a/requirements.txt b/requirements.txt
index e03d3c4f..0c0134a4 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,4 +2,5 @@ cmake
 cython>=0.28.0,<3.0.0
 numpy>=1.20.0
 ninja
+pytest
 scikit-build