diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py
index e9688423e..1ac8f353a 100644
--- a/elastica/rod/knot_theory.py
+++ b/elastica/rod/knot_theory.py
@@ -251,16 +251,29 @@ def _compute_twist(center_line, normal_collection):
         alpha /= 2 * np.pi
         gamma /= 2 * np.pi
         twist_temp = alpha + gamma
+
+        # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other.
+        idx = np.where(np.isnan(twist_temp))[0]
+        for i in idx:
+            # compute angle between two vectors
+            dot_product = np.dot(
+                projection_of_normal_collection[:, i],
+                projection_of_normal_collection[:, i + 1],
+            )
+            angle = np.arccos(dot_product)
+            cross = np.cross(
+                projection_of_normal_collection[:, i],
+                projection_of_normal_collection[:, i + 1],
+            )
+            angle *= np.sign(np.dot(tangent[:, i], cross))
+            twist_temp[i] = angle / (2 * np.pi)
+
         # Make sure twist is between (-1/2 to 1/2) as defined in pg 313 Klenin & Langowski 2000
         idx = np.where(twist_temp > 0.5)[0]
         twist_temp[idx] -= 1
         idx = np.where(twist_temp < -0.5)[0]
         twist_temp[idx] += 1
 
-        # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other.
-        idx = np.where(np.isnan(twist_temp))[0]
-        twist_temp[idx] = 0.0
-
         local_twist[k, :] = twist_temp
         total_twist[k] = np.sum(twist_temp)
 
diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py
index a5630d73a..6911fe0a2 100644
--- a/tests/test_rod/test_knot_theory.py
+++ b/tests/test_rod/test_knot_theory.py
@@ -265,3 +265,50 @@ def test_knot_theory_compute_additional_segment_none_case(
     assert_allclose(new_center_line.shape, [timesteps, 3, n_elem])
     assert_allclose(beginning_direction.shape[0], timesteps)
     assert_allclose(end_direction.shape[0], timesteps)
+
+
+@pytest.mark.parametrize("n_elem", [1, 2, 10, 50, 100])
+@pytest.mark.parametrize(
+    "diff_angle",
+    [
+        0.0,
+        np.pi * 0.1,
+        np.pi * 0.5,
+        np.pi * 0.9,
+        np.pi * 0.999,
+        -np.pi * 0.1,
+        -np.pi * 0.5,
+        -np.pi * 0.9,
+        -np.pi * 0.999,
+    ],
+)
+def test_knot_theory_compute_pure_twist(n_elem, diff_angle):
+    angle = np.arange(n_elem) * diff_angle
+
+    # setting up test params
+    normal = np.array([0.0, 1.0, 0.0])
+    base_length = 1.0
+
+    position_collection = np.zeros((3, n_elem + 1))
+    position_collection[2, :] = np.linspace(0, base_length, n_elem + 1)
+    normal_collection = np.zeros((3, n_elem))
+
+    for i in range(n_elem):
+        alpha = angle[i]
+        rot_matrix = np.array(
+            [
+                [np.cos(alpha), -np.sin(alpha), 0.0],
+                [np.sin(alpha), np.cos(alpha), 0.0],
+                [0.0, 0.0, 1.0],
+            ]
+        )
+        normal_temp = rot_matrix @ normal
+        normal_temp /= np.linalg.norm(normal_temp)
+
+        normal_collection[:, i] = normal_temp
+
+    total_twist, local_twist = compute_twist(
+        position_collection[None, ...], normal_collection[None, ...]
+    )
+    np.testing.assert_allclose(local_twist[0], np.diff(angle / (2 * np.pi)))
+    np.testing.assert_allclose(total_twist[0], angle[-1] / (2 * np.pi))