diff --git a/python/damask/_crystal.py b/python/damask/_crystal.py
index 04a276248..f5bcc2269 100644
--- a/python/damask/_crystal.py
+++ b/python/damask/_crystal.py
@@ -1241,4 +1241,4 @@ def relation_operations(self,
         o_p = np.stack((o.to_frame(uvw=o_p[:,0] if o_l != 'hP' else util.Bravais_to_Miller(uvtw=o_p[:,0])),
                         o.to_frame(hkl=o_p[:,1] if o_l != 'hP' else util.Bravais_to_Miller(hkil=o_p[:,1]))),
                         axis=-2)
-        return (o_l,Rotation.from_parallel(a=m_p,b=o_p))
+        return (o_l,Rotation.from_parallel(a=m_p,b=o_p,active=True))
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 3739ae936..026344948 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -811,7 +811,7 @@ def as_homochoric(self) -> np.ndarray:
         Returns
         -------
         h : numpy.ndarray, shape (...,3)
-            Homochoric vector (h_1, h_2, h_3) with ǀhǀ < (3/4*π)^(1/3).
+            Homochoric vector (h_1, h_2, h_3) with ǀhǀ < (3π/4)^(1/3).
 
         Examples
         --------
@@ -831,7 +831,7 @@ def as_cubochoric(self) -> np.ndarray:
         Returns
         -------
         x : numpy.ndarray, shape (...,3)
-            Cubochoric vector (x_1, x_2, x_3) with max(x_i) < 1/2*π^(2/3).
+            Cubochoric vector (x_1, x_2, x_3) with max(x_i) < 1/2 π^(2/3).
 
         Examples
         --------
@@ -1062,7 +1062,8 @@ def from_matrix(R: np.ndarray,
 
     @staticmethod
     def from_parallel(a: np.ndarray,
-                      b: np.ndarray ) -> 'Rotation':
+                      b: np.ndarray,
+                      active: bool = False ) -> 'Rotation':
         """
         Initialize from pairs of two orthogonal basis vectors.
 
@@ -1072,11 +1073,20 @@ def from_parallel(a: np.ndarray,
             Two three-dimensional vectors of first orthogonal basis.
         b : numpy.ndarray, shape (...,2,3)
             Corresponding three-dimensional vectors of second basis.
+        active : bool, optional
+            Consider rotations as active, i.e. return (B^-1⋅A) instead of (B⋅A^-1).
+            Defaults to False.
 
         Returns
         -------
         new : damask.Rotation
 
+        Notes
+        -----
+        If rotations $A = [a_1,a_2,a_1 × a_2]^T$ and B = $[b_1,b_2,b_1 × b_2]^T$
+        are considered "active", the resulting rotation will be $B^{-1}⋅A$ instead
+        of the default result $B⋅A^{-1}$.
+
         Examples
         --------
         >>> import damask
@@ -1095,10 +1105,10 @@ def from_parallel(a: np.ndarray,
 
         am = np.stack([          a_[...,0,:],
                                              a_[...,1,:],
-                        np.cross(a_[...,0,:],a_[...,1,:]) ],axis=-1)
+                        np.cross(a_[...,0,:],a_[...,1,:]) ],axis=-1 if active else -2)
         bm = np.stack([          b_[...,0,:],
                                              b_[...,1,:],
-                        np.cross(b_[...,0,:],b_[...,1,:]) ],axis=-1)
+                        np.cross(b_[...,0,:],b_[...,1,:]) ],axis=-1 if active else -2)
 
         return Rotation.from_basis(am).misorientation(Rotation.from_basis(bm))
 
@@ -1156,7 +1166,7 @@ def from_homochoric(h: np.ndarray,
         Parameters
         ----------
         h : numpy.ndarray, shape (...,3)
-            Homochoric vector (h_1, h_2, h_3) with ǀhǀ < (3/4*π)^(1/3).
+            Homochoric vector (h_1, h_2, h_3) with ǀhǀ < (3π/4)^(1/3).
         P : int ∈ {-1,1}, optional
             Sign convention. Defaults to -1.
 
@@ -1186,7 +1196,7 @@ def from_cubochoric(x: np.ndarray,
         Parameters
         ----------
         x : numpy.ndarray, shape (...,3)
-            Cubochoric vector (x_1, x_2, x_3) with max(x_i) < 1/2*π^(2/3).
+            Cubochoric vector (x_1, x_2, x_3) with max(x_i) < 1/2 π^(2/3).
         P : int ∈ {-1,1}, optional
             Sign convention. Defaults to -1.
 
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 4cb198741..8c4aa9b41 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -789,8 +789,16 @@ def test_matrix(self,multidim_rotations,normalize):
     def test_parallel(self,multidim_rotations):
         a = np.array([[1.42,0.0,0.0],
                       [0.0,0.3,0.0]])
-        assert Rotation.from_parallel(a,multidim_rotations[...,np.newaxis]@a).allclose( multidim_rotations)
-        assert Rotation.from_parallel(multidim_rotations[...,np.newaxis]@a,a).allclose(~multidim_rotations)
+        b = ~multidim_rotations[...,np.newaxis]@a            # actively rotate a axes as new b axes
+        assert Rotation.from_parallel(a,b).allclose( multidim_rotations)
+        assert Rotation.from_parallel(b,a).allclose(~multidim_rotations)
+
+
+    def test_basis_parallel_consistency(self,multidim_rotations):
+        M = multidim_rotations.as_matrix()
+        assert Rotation.from_basis(M).allclose(
+               Rotation.from_parallel(np.broadcast_to(np.identity(3),M.shape)[...,:2,:],
+                                      M[...,:2,:]))
 
 
     @pytest.mark.parametrize('normalize',[True,False])