diff --git a/orix/quaternion/quaternion.py b/orix/quaternion/quaternion.py
index 8b69a7f2..b89df394 100644
--- a/orix/quaternion/quaternion.py
+++ b/orix/quaternion/quaternion.py
@@ -222,7 +222,7 @@ def __mul__(self, other: Quaternion | Vector3d) -> Quaternion | Vector3d:
                 )
             else:  # pragma: no cover
                 v = qu_rotate_vec(self.unit.data, other.data)
-            if isinstance(other, Miller):
+            if hasattr(object, 'phase'):
                 m = other.__class__(xyz=v, phase=other.phase)
                 m.coordinate_format = other.coordinate_format
                 return m
diff --git a/orix/vector/miller.py b/orix/vector/miller.py
index b6a59c54..381b403e 100644
--- a/orix/vector/miller.py
+++ b/orix/vector/miller.py
@@ -21,6 +21,7 @@
 
 from copy import deepcopy
 from itertools import product
+import logging
 from typing import TYPE_CHECKING
 
 from diffpy.structure import Lattice
@@ -34,6 +35,9 @@
     from orix.quaternion.symmetry import Symmetry
 
 
+_logger = logging.getLogger(__name__)
+
+
 class Miller(Vector3d):
     r"""Direct crystal lattice vectors (uvw or UVTW) and reciprocal
     crystal lattice vectors (hkl or hkil), the latter known as Miller
@@ -366,6 +370,36 @@ def __getitem__(self, key) -> Miller:
         m.coordinate_format = self.coordinate_format
         return m
 
+    def __neg__(self) -> Vector3d:
+        out = self.__class__(-self.data, phase=self.phase)
+        out.coordinate_format = self.coordinate_format
+        return out
+
+    def __add__(
+        self, other: int | float | list | tuple | np.ndarray | Vector3d | Miller
+    ) -> Miller:
+        data = None
+        if isinstance(other, Vector3d):
+            self._check_for_crystal_mismatch(other)
+            if not  hasattr(other, "phase"):
+                # assume the other vector is being added in the same
+                # coordinate frame as the crystal vector
+                other = self.__class__(**{self.coordinate_format: other.data}, phase=self.phase)
+            data = self.data + other.data
+        elif isinstance(other, (int, float)):
+            data = self.data + other
+        elif isinstance(other, (list, tuple)):
+            other = np.array(other)
+        if isinstance(other, np.ndarray):
+            data = self.data - other[..., np.newaxis]
+
+        if data is None:
+            return NotImplemented
+
+        out = self.__class__(xyz=data, phase = self.phase)
+        out.coordinate_format = self.coordinate_format
+        return out
+
     # ------------------------ Class methods ------------------------- #
 
     @classmethod
diff --git a/orix/vector/vector3d.py b/orix/vector/vector3d.py
index c4d4caa3..d354ad89 100644
--- a/orix/vector/vector3d.py
+++ b/orix/vector/vector3d.py
@@ -20,6 +20,7 @@
 from __future__ import annotations
 
 from copy import deepcopy
+import logging
 from typing import TYPE_CHECKING, Any
 
 import dask.array as da
@@ -35,6 +36,8 @@
 if TYPE_CHECKING:  # pragma: no cover
     from orix.quaternion.symmetry import Symmetry
 
+_logger = logging.getLogger(__name__)
+
 
 class Vector3d(Object3d):
     r"""Three-dimensional vectors.
@@ -246,10 +249,18 @@ def __add__(
         self, other: int | float | list | tuple | np.ndarray | Vector3d
     ) -> Vector3d:
         if isinstance(other, Vector3d):
-            return self.__class__(self.data + other.data)
+            if hasattr(other, "phase"):
+                # If 'other' has a phase attribute (e.g, 'Miller'),
+                # use its addition function instead.
+                return other + self
+            else:
+                # Check for symmetry and reference fram mismatches
+                self._check_for_crystal_mismatch(other)
+                return self.__class__(self.data + other.data)
         elif isinstance(other, (int, float)):
             return self.__class__(self.data + other)
         elif isinstance(other, (list, tuple)):
+            # convert lists and tuples into numpy arrays
             other = np.array(other)
 
         if isinstance(other, np.ndarray):
@@ -257,41 +268,26 @@ def __add__(
 
         return NotImplemented
 
-    def __radd__(self, other: int | float | list | tuple | np.ndarray) -> Vector3d:
-        if isinstance(other, (int, float)):
-            return self.__class__(other + self.data)
-        elif isinstance(other, (list, tuple)):
-            other = np.array(other)
-
-        if isinstance(other, np.ndarray):
-            return self.__class__(other[..., np.newaxis] + self.data)
-
+    def __radd__(
+        self, other: int | float | list | tuple | np.ndarray
+    ) -> Vector3d:
+        if isinstance(other, (int, float, list, tuple, np.ndarray, Vector3d)):
+            return self + other
         return NotImplemented
 
     def __sub__(
         self, other: int | float | list | tuple | np.ndarray | Vector3d
     ) -> Vector3d:
-        if isinstance(other, Vector3d):
-            return self.__class__(self.data - other.data)
-        elif isinstance(other, (int, float)):
-            return self.__class__(self.data - other)
-        elif isinstance(other, (list, tuple)):
+        # Convert to `self + (-other)` to avoid code repetition
+        if isinstance(other, (list, tuple)):
             other = np.array(other)
+        return self + (-other)
 
-        if isinstance(other, np.ndarray):
-            return self.__class__(self.data - other[..., np.newaxis])
-
-        return NotImplemented
-
-    def __rsub__(self, other: int | float | list | tuple | np.ndarray) -> Vector3d:
-        if isinstance(other, (int, float)):
-            return self.__class__(other - self.data)
-        elif isinstance(other, (list, tuple)):
-            other = np.array(other)
-
-        if isinstance(other, np.ndarray):
-            return self.__class__(other[..., np.newaxis] - self.data)
-
+    def __rsub__(
+        self, other: int | float | list | tuple | np.ndarray
+    ) -> Vector3d:
+        if isinstance(other, (int, float, list, tuple, np.ndarray, Vector3d)):
+            return -self + other
         return NotImplemented
 
     def __mul__(
@@ -312,15 +308,11 @@ def __mul__(
 
         return NotImplemented
 
-    def __rmul__(self, other: int | float | list | tuple | np.ndarray) -> Vector3d:
-        if isinstance(other, (int, float)):
-            return self.__class__(other * self.data)
-        elif isinstance(other, (list, tuple)):
-            other = np.array(other)
-
-        if isinstance(other, np.ndarray):
-            return self.__class__(other[..., np.newaxis] * self.data)
-
+    def __rmul__(
+        self, other: int | float | list | tuple | np.ndarray
+    ) -> Vector3d:
+        if isinstance(other, (int, float, list, tuple, np.ndarray, Vector3d)):
+            return self * other
         return NotImplemented
 
     def __truediv__(
@@ -343,6 +335,7 @@ def __rtruediv__(self, other: Any) -> None:
 
     def __eq__(self, other: Any) -> np.ndarray:
         if isinstance(other, Vector3d):
+            self._check_for_symmetry_mismatch(other)
             return np.all(self.data == other.data, axis=-1)
         else:
             return self.data == other
@@ -1726,3 +1719,33 @@ def _dot_outer_dask(self, other: Vector3d, chunk_size: int = 20) -> da.Array:
         v2 = da.from_array(other.data, chunks=chunks2)
 
         return da.tensordot(v1, v2, axes=(v1.ndim - 1, v2.ndim - 1))
+
+    def _check_for_crystal_mismatch(self, other) -> None:
+        # NOTE: Vector3D does not include phase or coordinate_format
+        # information, but some of it's subclasses (most importantly, Miller)
+        # do. Adding the checker function to Vector3D allows logging
+        # consistent warnings regardless of operation order (ie, v+m or m+v).
+        phase_checks = [hasattr(self, "phase"), hasattr(other, "phase")]
+        if np.all(phase_checks):
+            if self.phase.point_group != other.phase.point_group:
+                pg1 = self.phase.point_group.name
+                pg2 = other.phase.point_group.name
+                _logger.warning(
+                    f"WARNING: Point group mismatch detected between crystal vectors ({pg1} versus {pg2})."
+                )
+        elif np.any(phase_checks):
+            _logger.warning(
+                "WARNING: Reference frame mismatch detected (Euclidean versus crystal vectors)."
+            )
+        crd_checks = [
+            hasattr(self, "coordinate_format"),
+            hasattr(other, "coordinate_format"),
+        ]
+        if np.all(crd_checks):
+            cf1 = self.coordinate_format
+            cf2 = other.coordinate_format
+            if cf1 != cf2:
+                _logger.warning(
+                    f"WARNING: Coordinates mismatch detected between crystal vectors ({cf1} versus {cf2})."
+                )
+        return None