diff --git a/ChangeLog.md b/ChangeLog.md
index f2b01ab44..7a3ef9fe8 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -15,6 +15,7 @@ and this project adheres to
 * The ``neighbors`` argument to ``env_neighbors`` for ``EnvironmentMotifMatch`` class.
 * The ``neighbors`` argument to ``cluster_neighbors`` for ``EnvironmentCluster`` class.
 * `freud.order.Nematic` uses orientation vectors instead of quaternions and a nematic director.
+* `freud.order.Nematic` raises a warning when the zero is vector passed.
 * Remove zero-padding from arrays in `freud.environment.EnvironmentCluster` and `freud.environment.EnvironmentMotifMatch` and replace with ragged lists of NumPy arrays.
 
 ### Removed
diff --git a/doc/source/reference/credits.rst b/doc/source/reference/credits.rst
index eb62b6265..5a09883a1 100644
--- a/doc/source/reference/credits.rst
+++ b/doc/source/reference/credits.rst
@@ -388,6 +388,7 @@ Kody Takada
 Alain Kadar
 
 * Introduced mass dependence for ``ClusterProperties`` class, inertia tensor and radius of gyration.
+* Added check for passing zero vector to ``Nematic`` class.
 
 Melody Zhang
 
diff --git a/freud/order.pyx b/freud/order.pyx
index 5a3212797..15e1c9162 100644
--- a/freud/order.pyx
+++ b/freud/order.pyx
@@ -208,11 +208,15 @@ cdef class Nematic(_Compute):
                 Orientation vectors for which to calculate the order parameter.
         """   # noqa: E501
         if orientations.shape[1] == 4:
-            raise ValueError('In freud versions >=3.0.0, Nematic.compute() takes'
+            raise ValueError('In freud versions >=3.0.0, Nematic.compute() takes '
                              '3d orientation vectors instead of 4d quaternions.')
         orientations = freud.util._convert_array(
             orientations, shape=(None, 3))
 
+        if len(np.where(~orientations.any(axis=1))[0])!=0:
+            warnings.warn('Including zero vector in the orientations array '
+                          'may lead to undefined behavior.',
+                          UserWarning)
         cdef const float[:, ::1] l_orientations = orientations
         cdef unsigned int num_particles = l_orientations.shape[0]
 
diff --git a/tests/test_order_Nematic.py b/tests/test_order_Nematic.py
index 6acbe1058..f462ac793 100644
--- a/tests/test_order_Nematic.py
+++ b/tests/test_order_Nematic.py
@@ -94,6 +94,25 @@ def test_imperfect(self):
         npt.assert_allclose(op_perp.nematic_tensor, np.diag([-0.5, 1, -0.5]), atol=1e-1)
         assert not np.all(op_perp.nematic_tensor == np.diag([-0.5, 1, -0.5]))
 
+    def test_warning(self):
+        """Test that supplying a zero orientation vector raises a warning."""
+        N = 10000
+        np.random.seed(0)
+        u = [1, 0, 0]
+
+        # Generate orientations close to the u
+        orientations = np.random.normal(
+            np.repeat(np.expand_dims(u, axis=0), repeats=N, axis=0), 0.1
+        )
+
+        # Change first orientation to zero vector
+        orientations[0] = np.array([0, 0, 0])
+
+        op = freud.order.Nematic()
+
+        with pytest.warns(UserWarning):
+            op.compute(orientations)
+
     def test_repr(self):
         op = freud.order.Nematic()
         assert str(op) == str(eval(repr(op)))