diff --git a/tests/conftest.py b/tests/conftest.py
index 8362e35a..b454edbb 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,12 +4,21 @@
 
 import gemmi
 import numpy as np
-import pandas as pd
 import pytest
 
 import reciprocalspaceship as rs
 
 
+@pytest.fixture
+def hkls():
+    """
+    Return all Miller indices with H, K, L values between [-5, 5]
+    """
+    hmin, hmax = -5, 5
+    H = np.mgrid[hmin : hmax + 1, hmin : hmax + 1, hmin : hmax + 1].reshape((3, -1)).T
+    return H
+
+
 @pytest.fixture
 def dataset_hkl():
     """
diff --git a/tests/utils/test_symop.py b/tests/utils/test_symop.py
new file mode 100644
index 00000000..d16d7329
--- /dev/null
+++ b/tests/utils/test_symop.py
@@ -0,0 +1,32 @@
+import gemmi
+import numpy as np
+import pytest
+
+import reciprocalspaceship as rs
+
+
+@pytest.mark.parametrize("H_even", [True, False])
+@pytest.mark.parametrize(
+    "op_str", ["x,y,z", "2*x,2*y,2*z", "x,z,y", "1/2*x,y,z", "1/2*x,1/2*y,1/2*z"]
+)
+def test_apply_to_hkl(hkls, H_even, op_str):
+    """
+    Test rs.utils.apply_to_hkl() detects symops that yield fractional Miller indices.
+
+    apply_to_hkl() should raise a RuntimeError if the combination of `H` and `op`
+    yield fractional Miller indices, and should return new Miller indices all other
+    cases.
+    """
+    if H_even:
+        hkls = hkls[~np.any(hkls % 2, axis=1)]
+
+    op = gemmi.Op(op_str)
+
+    if ((np.array(op.rot) / op.DEN) % 1 == 0).all() or H_even:
+        H_result = rs.utils.apply_to_hkl(hkls, op)
+        H_expected = np.array([op.apply_to_hkl(hkl) for hkl in hkls])
+        assert np.array_equal(H_expected, H_result)
+        assert H_result.dtype is np.dtype(np.int32)
+    else:
+        with pytest.raises(RuntimeError):
+            H_result = rs.utils.apply_to_hkl(hkls, op)