diff --git a/VERSION b/VERSION
index b7128c9af..c6d828650 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v3.0.0-alpha2-35-g1ebd10745
+v3.0.0-alpha2-153-gf8dd5df0c
diff --git a/python/damask/_config.py b/python/damask/_config.py
index 24245f4bd..76955588f 100644
--- a/python/damask/_config.py
+++ b/python/damask/_config.py
@@ -21,6 +21,9 @@ class NiceDumper(yaml.SafeDumper):
         return self.represent_data(dict(data)) if isinstance(data, dict) and type(data) != dict else \
                super().represent_data(data)
 
+    def ignore_aliases(self, data):
+        """No references."""
+        return True
 
 class Config(dict):
     """YAML-based configuration."""
diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py
index 6de2283f4..b94e9897a 100644
--- a/python/damask/_configmaterial.py
+++ b/python/damask/_configmaterial.py
@@ -9,6 +9,16 @@ from . import Orientation
 class ConfigMaterial(Config):
     """Material configuration."""
 
+    _defaults = {'material': [],
+                 'homogenization': {},
+                 'phase': {}}
+
+    def __init__(self,d={}):
+        """Initialize object with default dictionary keys."""
+        super().__init__(d)
+        for k,v in self._defaults.items():
+            if k not in self: self[k] = v
+
     def save(self,fname='material.yaml',**kwargs):
         """
         Save to yaml file.
@@ -75,6 +85,8 @@ class ConfigMaterial(Config):
                 fraction: 1.0
                 phase: Steel
             homogenization: SX
+        homogenization: {}
+        phase: {}
 
         """
         constituents_ = {k:table.get(v) for k,v in constituents.items()}
@@ -261,6 +273,8 @@ class ConfigMaterial(Config):
                 fraction: 1.0
                 phase: Aluminum
             homogenization: SX
+        homogenization: {}
+        phase: {}
 
         """
         length = -1
@@ -274,6 +288,7 @@ class ConfigMaterial(Config):
 
         c = [{} for _ in range(length)] if constituents is None else \
             [{'constituents':u} for u in ConfigMaterial._constituents(**constituents)]
+
         if len(c) == 1: c = [copy.deepcopy(c[0]) for _ in range(length)]
 
         if length != 1 and length != len(c):
diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index 05301561f..e6434813e 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -226,9 +226,9 @@ class Orientation(Rotation):
 
         """
         return  super().__eq__(other) \
-            and self.family     == other.family \
-            and self.lattice    == other.lattice \
-            and self.parameters == other.parameters
+            and hasattr(other, 'family')     and self.family     == other.family \
+            and hasattr(other, 'lattice')    and self.lattice    == other.lattice \
+            and hasattr(other, 'parameters') and self.parameters == other.parameters
 
 
     def __matmul__(self,other):
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index f4d2cf248..78029b59a 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -144,6 +144,11 @@ class Rotation:
         return self.copy(rotation=Rotation(np.block([np.cos(pwr*phi),np.sin(pwr*phi)*p]))._standardize())
 
 
+    def __mul__(self,other):
+        """Standard multiplication is not implemented."""
+        raise NotImplementedError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
+
+
     def __matmul__(self,other):
         """
         Rotation of vector, second or fourth order tensor, or rotation object.
@@ -199,8 +204,16 @@ class Rotation:
 
 
     def append(self,other):
-        """Extend rotation array along first dimension with other array."""
-        return self.copy(rotation=np.vstack((self.quaternion,other.quaternion)))
+        """
+        Extend rotation array along first dimension with other array(s).
+
+        Parameters
+        ----------
+            other : Rotation or list of Rotations.
+
+        """
+        return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion,
+                                                      [self]+other if isinstance(other,list) else [self,other]))))
 
 
     def flatten(self,order = 'C'):
@@ -258,7 +271,7 @@ class Rotation:
             """Intermediate representation supporting quaternion averaging."""
             return np.einsum('...i,...j',quat,quat)
 
-        if not weights:
+        if weights is None:
             weights = np.ones(self.shape,dtype=float)
 
         eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights[...,np.newaxis,np.newaxis],axis=-3) \
diff --git a/python/tests/reference/ConfigMaterial/material.yaml b/python/tests/reference/ConfigMaterial/material.yaml
index 933e295b3..fbba6a631 100644
--- a/python/tests/reference/ConfigMaterial/material.yaml
+++ b/python/tests/reference/ConfigMaterial/material.yaml
@@ -1,10 +1,10 @@
 homogenization:
   SX:
-    N_constituents: 2
-    mech: {type: none}
+    N_constituents: 1
+    mechanics: {type: none}
   Taylor:
     N_constituents: 2
-    mech: {type: isostrain}
+    mechanics: {type: isostrain}
 
 material:
   - constituents:
@@ -34,11 +34,11 @@ material:
 phase:
   Aluminum:
     lattice: cF
-    mech:
+    mechanics:
       output: [F, P, F_e, F_p, L_p]
       elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
   Steel:
     lattice: cI
-    mech:
+    mechanics:
       output: [F, P, F_e, F_p, L_p]
       elasticity: {C_11: 233.3e9, C_12: 135.5e9, C_44: 118.0e9, type: hooke}
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index c60029046..36e3a3ac9 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -800,6 +800,14 @@ class TestRotation:
         print(f'append 2x {shape} --> {s.shape}')
         assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
 
+    @pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)])
+    def test_append_list(self,shape):
+        r = Rotation.from_random(shape=shape)
+        p = Rotation.from_random(shape=shape)
+        s = r.append([r,p])
+        print(f'append 3x {shape} --> {s.shape}')
+        assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
+
     @pytest.mark.parametrize('quat,standardized',[
                                                   ([-1,0,0,0],[1,0,0,0]),
                                                   ([-0.5,-0.5,-0.5,-0.5],[0.5,0.5,0.5,0.5]),
diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90
index d8ab6390d..674ec4c43 100644
--- a/src/commercialFEM_fileList.f90
+++ b/src/commercialFEM_fileList.f90
@@ -11,7 +11,6 @@
 #include "config.f90"
 #include "LAPACK_interface.f90"
 #include "math.f90"
-#include "quaternions.f90"
 #include "rotations.f90"
 #include "element.f90"
 #include "HDF5_utilities.f90"
diff --git a/src/constitutive.f90 b/src/constitutive.f90
index 667a23127..7dac3e129 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -859,10 +859,11 @@ subroutine crystallite_init
     co, &                                                                                           !< counter in integration point component loop
     ip, &                                                                                           !< counter in integration point loop
     el, &                                                                                           !< counter in element loop
+    so, &
     cMax, &                                                                                         !< maximum number of  integration point components
     iMax, &                                                                                         !< maximum number of integration points
-    eMax, &                                                                                         !< maximum number of elements
-    so
+    eMax                                                                                            !< maximum number of elements
+   
 
   class(tNode), pointer :: &
     num_crystallite, &
diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90
index f5a5cd0a2..44e777c80 100644
--- a/src/constitutive_mech.f90
+++ b/src/constitutive_mech.f90
@@ -1109,6 +1109,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
 
   enddo iteration
 
+
   contains
 
   !--------------------------------------------------------------------------------------------------
diff --git a/src/quaternions.f90 b/src/quaternions.f90
deleted file mode 100644
index c5c43e3c1..000000000
--- a/src/quaternions.f90
+++ /dev/null
@@ -1,534 +0,0 @@
-!---------------------------------------------------------------------------------------------------
-!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Philip Eisenlohr, Michigan State University
-!> @brief general quaternion math, not limited to unit quaternions
-!> @details w is the real part, (x, y, z) are the imaginary parts.
-!> @details https://en.wikipedia.org/wiki/Quaternion
-!---------------------------------------------------------------------------------------------------
-module quaternions
-  use prec
-
-  implicit none
-  private
-
-  real(pReal), parameter, public :: P = -1.0_pReal                                                  !< parameter for orientation conversion.
-
-  type, public :: quaternion
-    real(pReal), private :: w = 0.0_pReal
-    real(pReal), private :: x = 0.0_pReal
-    real(pReal), private :: y = 0.0_pReal
-    real(pReal), private :: z = 0.0_pReal
-
-
-  contains
-    procedure, private :: add__
-    procedure, private :: pos__
-    generic,   public  :: operator(+) => add__,pos__
-
-    procedure, private :: sub__
-    procedure, private :: neg__
-    generic,   public  :: operator(-) => sub__,neg__
-
-    procedure, private :: mul_quat__
-    procedure, private :: mul_scal__
-    generic,   public  :: operator(*) => mul_quat__, mul_scal__
-
-    procedure, private :: div_quat__
-    procedure, private :: div_scal__
-    generic,   public  :: operator(/) => div_quat__, div_scal__
-
-    procedure, private :: eq__
-    generic,   public  :: operator(==) => eq__
-
-    procedure, private :: neq__
-    generic,   public  :: operator(/=) => neq__
-
-    procedure, private :: pow_quat__
-    procedure, private :: pow_scal__
-    generic,   public  :: operator(**) => pow_quat__, pow_scal__
-
-    procedure, public  :: abs   => abs__
-    procedure, public  :: conjg => conjg__
-    procedure, public  :: real  => real__
-    procedure, public  :: aimag => aimag__
-
-    procedure, public  :: homomorphed
-    procedure, public  :: asArray
-    procedure, public  :: inverse
-
-  end type
-
-  interface assignment (=)
-    module procedure assign_quat__
-    module procedure assign_vec__
-  end interface assignment (=)
-
-  interface quaternion
-    module procedure init__
-  end interface quaternion
-
-  interface abs
-    procedure abs__
-  end interface abs
-
-  interface dot_product
-    procedure dot_product__
-  end interface dot_product
-
-  interface conjg
-    module procedure conjg__
-  end interface conjg
-
-  interface exp
-    module procedure exp__
-  end interface exp
-
-  interface log
-    module procedure log__
-  end interface log
-
-  interface real
-    module procedure real__
-  end interface real
-
-  interface aimag
-    module procedure aimag__
-  end interface aimag
-
-  public :: &
-    quaternions_init, &
-    assignment(=), &
-    conjg, aimag, &
-    log, exp, &
-    abs, dot_product, &
-    inverse, &
-    real
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Do self test.
-!--------------------------------------------------------------------------------------------------
-subroutine quaternions_init
-
-  print'(/,a)', ' <<<+-  quaternions init  -+>>>'; flush(6)
-
-  call selfTest
-
-end subroutine quaternions_init
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief construct a quaternion from a 4-vector
-!---------------------------------------------------------------------------------------------------
-type(quaternion) pure function init__(array)
-
-  real(pReal), intent(in), dimension(4) :: array
-
-  init__%w = array(1)
-  init__%x = array(2)
-  init__%y = array(3)
-  init__%z = array(4)
-
-end function init__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief assign a quaternion
-!---------------------------------------------------------------------------------------------------
-elemental pure subroutine assign_quat__(self,other)
-
-  type(quaternion), intent(out) :: self
-  type(quaternion), intent(in)  :: other
-
-  self = [other%w,other%x,other%y,other%z]
-
-end subroutine assign_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief assign a 4-vector
-!---------------------------------------------------------------------------------------------------
-pure subroutine assign_vec__(self,other)
-
-  type(quaternion), intent(out)                :: self
-  real(pReal),       intent(in), dimension(4)  :: other
-
-  self%w = other(1)
-  self%x = other(2)
-  self%y = other(3)
-  self%z = other(4)
-
-end subroutine assign_vec__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief add a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function add__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  add__ = [ self%w,  self%x,  self%y ,self%z] &
-        + [other%w, other%x, other%y,other%z]
-
-end function add__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief return (unary positive operator)
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pos__(self)
-
-  class(quaternion), intent(in) :: self
-
-  pos__ = self * (+1.0_pReal)
-
-end function pos__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief subtract a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function sub__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  sub__ = [ self%w,  self%x,  self%y ,self%z] &
-        - [other%w, other%x, other%y,other%z]
-
-end function sub__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief negate (unary negative operator)
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function neg__(self)
-
-  class(quaternion), intent(in) :: self
-
-  neg__ = self * (-1.0_pReal)
-
-end function neg__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief multiply with a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function mul_quat__(self,other)
-
-  class(quaternion), intent(in) :: self, other
-
-  mul_quat__%w = self%w*other%w - self%x*other%x -      self%y*other%y - self%z*other%z
-  mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y)
-  mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z)
-  mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x)
-
-end function mul_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief multiply with a scalar
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function mul_scal__(self,scal)
-
-  class(quaternion), intent(in) :: self
-  real(pReal),       intent(in) :: scal
-
-  mul_scal__ = [self%w,self%x,self%y,self%z]*scal
-
-end function mul_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief divide by a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function div_quat__(self,other)
-
-  class(quaternion), intent(in) :: self, other
-
-  div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
-
-end function div_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief divide by a scalar
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function div_scal__(self,scal)
-
-  class(quaternion), intent(in) :: self
-  real(pReal),       intent(in) :: scal
-
-  div_scal__ = [self%w,self%x,self%y,self%z]/scal
-
-end function div_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief test equality
-!---------------------------------------------------------------------------------------------------
-logical elemental pure function eq__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  eq__ = all(dEq([ self%w, self%x, self%y, self%z], &
-                 [other%w,other%x,other%y,other%z]))
-
-end function eq__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief test inequality
-!---------------------------------------------------------------------------------------------------
-logical elemental pure function neq__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  neq__ = .not. self%eq__(other)
-
-end function neq__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief raise to the power of a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pow_quat__(self,expon)
-
-  class(quaternion), intent(in) :: self
-  type(quaternion),  intent(in) :: expon
-
-  pow_quat__ = exp(log(self)*expon)
-
-end function pow_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief raise to the power of a scalar
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pow_scal__(self,expon)
-
-  class(quaternion), intent(in) :: self
-  real(pReal),       intent(in) :: expon
-
-  pow_scal__ = exp(log(self)*expon)
-
-end function pow_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief take exponential
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function exp__(a)
-
-  class(quaternion), intent(in) :: a
-  real(pReal)                   :: absImag
-
-  absImag = norm2(aimag(a))
-
-  exp__ = merge(exp(a%w) * [               cos(absImag), &
-                             a%x/absImag * sin(absImag), &
-                             a%y/absImag * sin(absImag), &
-                             a%z/absImag * sin(absImag)], &
-                IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
-                dNeq0(absImag))
-
-end function exp__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief take logarithm
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function log__(a)
-
-  class(quaternion), intent(in) :: a
-  real(pReal)                   :: absImag
-
-  absImag = norm2(aimag(a))
-
-  log__ = merge([log(abs(a)), &
-                 a%x/absImag * acos(a%w/abs(a)), &
-                 a%y/absImag * acos(a%w/abs(a)), &
-                 a%z/absImag * acos(a%w/abs(a))], &
-                IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
-                dNeq0(absImag))
-
-end function log__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief return norm
-!---------------------------------------------------------------------------------------------------
-real(pReal) elemental pure function abs__(self)
-
-  class(quaternion), intent(in) :: self
-
-  abs__ = norm2([self%w,self%x,self%y,self%z])
-
-end function abs__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief calculate dot product
-!---------------------------------------------------------------------------------------------------
-real(pReal) elemental pure function dot_product__(a,b)
-
-  class(quaternion), intent(in) :: a,b
-
-  dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
-
-end function dot_product__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief take conjugate complex
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function conjg__(self)
-
-  class(quaternion), intent(in) :: self
-
-  conjg__ = [self%w,-self%x,-self%y,-self%z]
-
-end function conjg__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief homomorph
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function homomorphed(self)
-
-  class(quaternion), intent(in) :: self
-
-  homomorphed = - self
-
-end function homomorphed
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief return as plain array
-!---------------------------------------------------------------------------------------------------
-pure function asArray(self)
-
-  real(pReal), dimension(4)     :: asArray
-  class(quaternion), intent(in) :: self
-
-  asArray = [self%w,self%x,self%y,self%z]
-
-end function asArray
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief real part (scalar)
-!---------------------------------------------------------------------------------------------------
-pure function real__(self)
-
-  real(pReal)                   :: real__
-  class(quaternion), intent(in) :: self
-
-  real__ = self%w
-
-end function real__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief imaginary part (3-vector)
-!---------------------------------------------------------------------------------------------------
-pure function aimag__(self)
-
-  real(pReal), dimension(3)     :: aimag__
-  class(quaternion), intent(in) :: self
-
-  aimag__ = [self%x,self%y,self%z]
-
-end function aimag__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief inverse
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function inverse(self)
-
-  class(quaternion), intent(in) :: self
-
-  inverse = conjg(self)/abs(self)**2.0_pReal
-
-end function inverse
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief check correctness of some quaternions functions
-!--------------------------------------------------------------------------------------------------
-subroutine selfTest
-
-  real(pReal), dimension(4) :: qu
-  type(quaternion)          :: q, q_2
-
-  if(dNeq(abs(P),1.0_pReal))                           error stop 'P not in {-1,+1}'
-
-  call random_number(qu)
-  qu = (qu-0.5_pReal) * 2.0_pReal
-  q  = quaternion(qu)
-
-  q_2= qu
-  if(any(dNeq(q%asArray(),q_2%asArray())))             error stop 'assign_vec__'
-
-  q_2 = q + q
-  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))            error stop 'add__'
-
-  q_2 = q - q
-  if(any(dNeq0(q_2%asArray())))                        error stop 'sub__'
-
-  q_2 = q * 5.0_pReal
-  if(any(dNeq(q_2%asArray(),5.0_pReal*qu)))            error stop 'mul__'
-
-  q_2 = q / 0.5_pReal
-  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))            error stop 'div__'
-
-  q_2 = q * 0.3_pReal
-  if(dNeq0(abs(q)) .and. q_2 == q)                     error stop 'eq__'
-
-  q_2 = q
-  if(q_2 /= q)                                         error stop 'neq__'
-
-  if(dNeq(abs(q),norm2(qu)))                           error stop 'abs__'
-  if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) &
-                                                       error stop 'abs__/*conjg'
-
-  if(any(dNeq(q%asArray(),qu)))                        error stop 'eq__'
-  if(dNeq(q%real(),       qu(1)))                      error stop 'real()'
-  if(any(dNeq(q%aimag(),  qu(2:4))))                   error stop 'aimag()'
-
-  q_2 = q%homomorphed()
-  if(q                 /= q_2*    (-1.0_pReal))        error stop 'homomorphed'
-  if(dNeq(q_2%real(),     qu(1)*  (-1.0_pReal)))       error stop 'homomorphed/real'
-  if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal))))      error stop 'homomorphed/aimag'
-
-  q_2 = conjg(q)
-  if(dNeq(abs(q),abs(q_2)))                            error stop 'conjg/abs'
-  if(q /= conjg(q_2))                                  error stop 'conjg/involution'
-  if(dNeq(q_2%real(),     q%real()))                   error stop 'conjg/real'
-  if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal))))    error stop 'conjg/aimag'
-
-  if(abs(q) > 0.0_pReal) then
-    q_2 = q * q%inverse()
-    if(     dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal))  error stop 'inverse/real'
-    if(any(dNeq0(aimag(q_2),          1.0e-15_pReal))) error stop 'inverse/aimag'
-
-    q_2 = q/abs(q)
-    q_2 = conjg(q_2) - inverse(q_2)
-    if(any(dNeq0(q_2%asArray(),1.0e-15_pReal)))        error stop 'inverse/conjg'
-  endif
-  if(dNeq(dot_product(qu,qu),dot_product(q,q)))        error stop 'dot_product'
-
-#if !(defined(__GFORTRAN__) &&  __GNUC__ < 9)
-  if (norm2(aimag(q)) > 0.0_pReal) then
-    if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal))       error stop 'exp/log'
-    if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal))       error stop 'log/exp'
-  endif
-#endif
-
-end subroutine selfTest
-
-
-end module quaternions
diff --git a/src/rotations.f90 b/src/rotations.f90
index 888e73762..57dd16b53 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -47,16 +47,16 @@
 !---------------------------------------------------------------------------------------------------
 
 module rotations
-  use prec
   use IO
   use math
-  use quaternions
 
   implicit none
   private
 
+  real(pReal), parameter :: P = -1.0_pReal                                                          !< parameter for orientation conversion.
+
   type, public :: rotation
-    type(quaternion) :: q
+    real(pReal), dimension(4) :: q
     contains
       procedure, public :: asQuaternion
       procedure, public :: asEulers
@@ -103,7 +103,6 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine rotations_init
 
-  call quaternions_init
   print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
@@ -122,7 +121,7 @@ pure function asQuaternion(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(4)   :: asQuaternion
 
-  asQuaternion = self%q%asArray()
+  asQuaternion = self%q
 
 end function asQuaternion
 !---------------------------------------------------------------------------------------------------
@@ -131,7 +130,7 @@ pure function asEulers(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(3)   :: asEulers
 
-  asEulers = qu2eu(self%q%asArray())
+  asEulers = qu2eu(self%q)
 
 end function asEulers
 !---------------------------------------------------------------------------------------------------
@@ -140,7 +139,7 @@ pure function asAxisAngle(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(4)   :: asAxisAngle
 
-  asAxisAngle = qu2ax(self%q%asArray())
+  asAxisAngle = qu2ax(self%q)
 
 end function asAxisAngle
 !---------------------------------------------------------------------------------------------------
@@ -149,7 +148,7 @@ pure function asMatrix(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(3,3) :: asMatrix
 
-  asMatrix = qu2om(self%q%asArray())
+  asMatrix = qu2om(self%q)
 
 end function asMatrix
 !---------------------------------------------------------------------------------------------------
@@ -158,7 +157,7 @@ pure function asRodrigues(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(4)   :: asRodrigues
 
-  asRodrigues = qu2ro(self%q%asArray())
+  asRodrigues = qu2ro(self%q)
 
 end function asRodrigues
 !---------------------------------------------------------------------------------------------------
@@ -167,7 +166,7 @@ pure function asHomochoric(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(3)   :: asHomochoric
 
-  asHomochoric = qu2ho(self%q%asArray())
+  asHomochoric = qu2ho(self%q)
 
 end function asHomochoric
 
@@ -259,7 +258,7 @@ pure elemental function rotRot__(self,R) result(rRot)
   type(rotation)              :: rRot
   class(rotation), intent(in) :: self,R
 
-  rRot = rotation(self%q*R%q)
+  rRot = rotation(multiply_quaternion(self%q,R%q))
   call rRot%standardize()
 
 end function rotRot__
@@ -272,14 +271,14 @@ pure elemental subroutine standardize(self)
 
   class(rotation), intent(inout) :: self
 
-  if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed()
+  if (self%q(1) < 0.0_pReal) self%q = - self%q
 
 end subroutine standardize
 
 
 !---------------------------------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
-!> @brief rotate a vector passively (default) or actively
+!> @brief Rotate a vector passively (default) or actively.
 !---------------------------------------------------------------------------------------------------
 pure function rotVector(self,v,active) result(vRot)
 
@@ -288,9 +287,8 @@ pure function rotVector(self,v,active) result(vRot)
   real(pReal),     intent(in), dimension(3) :: v
   logical,         intent(in), optional     :: active
 
-  real(pReal),    dimension(3) :: v_normed
-  type(quaternion)             :: q
-  logical                      :: passive
+  real(pReal), dimension(4) :: v_normed, q
+  logical                   :: passive
 
   if (present(active)) then
     passive = .not. active
@@ -301,13 +299,13 @@ pure function rotVector(self,v,active) result(vRot)
   if (dEq0(norm2(v))) then
     vRot = v
   else
-    v_normed = v/norm2(v)
+    v_normed = [0.0_pReal,v]/norm2(v)
     if (passive) then
-      q = self%q * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) )
+      q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q)))
     else
-      q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
+      q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q))
     endif
-    vRot = q%aimag()*norm2(v)
+    vRot = q(2:4)*norm2(v)
   endif
 
 end function rotVector
@@ -315,8 +313,8 @@ end function rotVector
 
 !---------------------------------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
-!> @brief rotate a rank-2 tensor passively (default) or actively
-!> @details: rotation is based on rotation matrix
+!> @brief Rotate a rank-2 tensor passively (default) or actively.
+!> @details: Rotation is based on rotation matrix
 !---------------------------------------------------------------------------------------------------
 pure function rotTensor2(self,T,active) result(tRot)
 
@@ -403,7 +401,7 @@ pure elemental function misorientation(self,other)
   type(rotation)              :: misorientation
   class(rotation), intent(in) :: self, other
 
-  misorientation%q = other%q * conjg(self%q)
+  misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
 
 end function misorientation
 
@@ -1338,7 +1336,7 @@ end function cu2ho
 !--------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief determine to which pyramid a point in a cubic grid belongs
+!> @brief Determine to which pyramid a point in a cubic grid belongs.
 !--------------------------------------------------------------------------
 pure function GetPyramidOrder(xyz)
 
@@ -1362,7 +1360,39 @@ end function GetPyramidOrder
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief check correctness of some rotations functions
+!> @brief Multiply two quaternions.
+!--------------------------------------------------------------------------------------------------
+pure function multiply_quaternion(qu1,qu2)
+
+  real(pReal), dimension(4), intent(in) :: qu1, qu2
+  real(pReal), dimension(4) :: multiply_quaternion
+
+
+  multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) -      qu1(3)*qu2(3) - qu1(4)*qu2(4)
+  multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3))
+  multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4))
+  multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2))
+
+end function multiply_quaternion
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculate conjugate complex of a quaternion.
+!--------------------------------------------------------------------------------------------------
+pure function conjugate_quaternion(qu)
+
+  real(pReal), dimension(4), intent(in) :: qu
+  real(pReal), dimension(4) :: conjugate_quaternion
+
+
+  conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)]
+
+
+end function conjugate_quaternion
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Check correctness of some rotations functions.
 !--------------------------------------------------------------------------------------------------
 subroutine selfTest
 
@@ -1374,7 +1404,8 @@ subroutine selfTest
   real    :: A,B
   integer :: i
 
-  do i=1,10
+
+  do i = 1, 10
 
 #if defined(__GFORTRAN__) && __GNUC__<9
     if(i<7) cycle