diff --git a/VERSION b/VERSION
index 1485aa8b8..346a2ab96 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-3.0.0-alpha6-285-g4d131ec7b
+3.0.0-alpha6-291-g7d18ed307
diff --git a/src/lattice.f90 b/src/lattice.f90
index 79e3c1615..afff2a309 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -1966,8 +1966,8 @@ end function buildCoordinateSystem
 !--------------------------------------------------------------------------------------------------
 !> @brief Helper function to define transformation systems
 ! Needed to calculate Schmid matrix and rotated stiffness matrices.
-! @details: use c/a for cF -> cI transformation
-!           use a_cX for cF -> hP transformation
+! @details: use c/a for cF -> hP transformation
+!           use a_cX for cF -> cI transformation
 !--------------------------------------------------------------------------------------------------
 subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
 
diff --git a/src/math.f90 b/src/math.f90
index 75b9ddea2..94ad2c48d 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -127,8 +127,10 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
 
   integer, dimension(:,:), intent(inout) :: a
   integer, intent(in),optional :: istart,iend, sortDim
+
   integer :: ipivot,s,e,d
 
+
   if (present(istart)) then
     s = istart
   else
@@ -164,9 +166,11 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
     integer, dimension(:,:), intent(inout) :: a
     integer,                 intent(out)   :: p                                                     ! Pivot element
     integer,                 intent(in)    :: istart,iend,sort
-    integer, dimension(size(a,1))          :: tmp
+
+    integer, dimension(size(a,1)) :: tmp
     integer :: i,j
 
+
     do
       ! find the first element on the right side less than or equal to the pivot point
       do j = iend, istart, -1
@@ -204,8 +208,10 @@ pure function math_expand(what,how)
   real(pReal),   dimension(:), intent(in) :: what
   integer,       dimension(:), intent(in) :: how
   real(pReal), dimension(sum(how)) ::  math_expand
+
   integer :: i
 
+
   if (sum(how) == 0) return
 
   do i = 1, size(how)
@@ -221,9 +227,11 @@ end function math_expand
 pure function math_range(N)
 
   integer, intent(in) :: N                                                                          !< length of range
-  integer :: i
   integer, dimension(N) :: math_range
 
+  integer :: i
+
+
   math_range = [(i,i=1,N)]
 
 end function math_range
@@ -235,9 +243,11 @@ end function math_range
 pure function math_eye(d)
 
   integer, intent(in) :: d                                                                          !< tensor dimension
-  integer :: i
   real(pReal), dimension(d,d) :: math_eye
 
+  integer :: i
+
+
   math_eye = 0.0_pReal
   do i=1,d
     math_eye(i,i) = 1.0_pReal
@@ -302,6 +312,7 @@ real(pReal) pure function math_delta(i,j)
 
   integer, intent (in) :: i,j
 
+
   math_delta = merge(0.0_pReal, 1.0_pReal, i /= j)
 
 end function math_delta
@@ -315,6 +326,7 @@ pure function math_cross(A,B)
   real(pReal), dimension(3), intent(in) ::  A,B
   real(pReal), dimension(3) :: math_cross
 
+
   math_cross = [ A(2)*B(3) -A(3)*B(2), &
                  A(3)*B(1) -A(1)*B(3), &
                  A(1)*B(2) -A(2)*B(1) ]
@@ -329,6 +341,7 @@ pure function math_outer(A,B)
 
   real(pReal), dimension(:), intent(in) ::  A,B
   real(pReal), dimension(size(A,1),size(B,1)) ::  math_outer
+
   integer :: i,j
 
 
@@ -351,6 +364,7 @@ real(pReal) pure function math_inner(A,B)
   real(pReal), dimension(:),         intent(in) :: A
   real(pReal), dimension(size(A,1)), intent(in) :: B
 
+
   math_inner = sum(A*B)
 
 end function math_inner
@@ -363,6 +377,7 @@ real(pReal) pure function math_tensordot(A,B)
 
   real(pReal), dimension(3,3), intent(in) :: A,B
 
+
   math_tensordot = sum(A*B)
 
 end function math_tensordot
@@ -376,6 +391,7 @@ pure function math_mul3333xx33(A,B)
   real(pReal), dimension(3,3,3,3), intent(in) :: A
   real(pReal), dimension(3,3),     intent(in) :: B
   real(pReal), dimension(3,3) :: math_mul3333xx33
+
   integer :: i,j
 
 
@@ -395,11 +411,12 @@ end function math_mul3333xx33
 !--------------------------------------------------------------------------------------------------
 pure function math_mul3333xx3333(A,B)
 
-  integer :: i,j,k,l
   real(pReal), dimension(3,3,3,3), intent(in) :: A
   real(pReal), dimension(3,3,3,3), intent(in) :: B
   real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
 
+  integer :: i,j,k,l
+
 
 #ifndef __INTEL_COMPILER
   do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
@@ -424,6 +441,7 @@ pure function math_exp33(A,n)
   real(pReal) :: invFac
   integer     :: n_,i
 
+
   if (present(n)) then
     n_ = n
   else
@@ -456,6 +474,7 @@ pure function math_inv33(A)
   real(pReal) :: DetA
   logical     :: error
 
+
   call math_invert33(math_inv33,DetA,error,A)
   if (error) math_inv33 = 0.0_pReal
 
@@ -474,6 +493,7 @@ pure subroutine math_invert33(InvA, DetA, error, A)
   logical,                     intent(out) :: error
   real(pReal), dimension(3,3), intent(in)  :: A
 
+
   InvA(1,1) =  A(2,2) * A(3,3) - A(2,3) * A(3,2)
   InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1)
   InvA(3,1) =  A(2,1) * A(3,2) - A(2,2) * A(3,1)
@@ -504,7 +524,7 @@ end subroutine math_invert33
 !--------------------------------------------------------------------------------------------------
 pure function math_invSym3333(A)
 
-  real(pReal),dimension(3,3,3,3)            :: math_invSym3333
+  real(pReal),dimension(3,3,3,3) :: math_invSym3333
 
   real(pReal),dimension(3,3,3,3),intent(in) :: A
 
@@ -513,6 +533,7 @@ pure function math_invSym3333(A)
   real(pReal), dimension(6*6) :: work
   integer                     :: ierr_i, ierr_f
 
+
   temp66 = math_sym3333to66(A)
   call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
   call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
@@ -538,6 +559,7 @@ pure subroutine math_invert(InvA, error, A)
   real(pReal), dimension(size(A,1)**2) :: work
   integer                              :: ierr
 
+
   invA = A
   call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
   error = (ierr /= 0)
@@ -555,6 +577,7 @@ pure function math_symmetric33(m)
   real(pReal), dimension(3,3) :: math_symmetric33
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_symmetric33 = 0.5_pReal * (m + transpose(m))
 
 end function math_symmetric33
@@ -568,6 +591,7 @@ pure function math_skew33(m)
   real(pReal), dimension(3,3) :: math_skew33
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_skew33 = m - math_symmetric33(m)
 
 end function math_skew33
@@ -581,6 +605,7 @@ pure function math_spherical33(m)
   real(pReal), dimension(3,3) :: math_spherical33
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_spherical33 = math_I3 * math_trace33(m)/3.0_pReal
 
 end function math_spherical33
@@ -594,6 +619,7 @@ pure function math_deviatoric33(m)
   real(pReal), dimension(3,3) :: math_deviatoric33
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_deviatoric33 = m - math_spherical33(m)
 
 end function math_deviatoric33
@@ -606,6 +632,7 @@ real(pReal) pure function math_trace33(m)
 
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_trace33 = m(1,1) + m(2,2) + m(3,3)
 
 end function math_trace33
@@ -618,6 +645,7 @@ real(pReal) pure function math_det33(m)
 
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
              - m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) &
              + m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1))
@@ -632,6 +660,7 @@ real(pReal) pure function math_detSym33(m)
 
   real(pReal), dimension(3,3), intent(in) :: m
 
+
   math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) &
                   + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3)
 
@@ -761,6 +790,7 @@ pure function math_99to3333(m99)
 
   integer :: i,j
 
+
 #ifndef __INTEL_COMPILER
   do concurrent(i=1:9, j=1:9)
     math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
@@ -1012,6 +1042,7 @@ pure subroutine math_eigh33(w,v,m)
   real(pReal) :: T, U, norm, threshold
   logical :: error
 
+
   w = math_eigvalsh33(m)
 
   v(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
@@ -1066,6 +1097,7 @@ pure function math_rotationalPart(F) result(R)
     I_F                                                                                             ! first two invariants of F
   real(pReal) :: x,Phi
 
+
   C = matmul(transpose(F),F)
   I_C = math_invariantsSym33(C)
   I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))]
@@ -1105,6 +1137,7 @@ pure function math_eigvalsh(m)
   integer :: ierr
   real(pReal), dimension(size(m,1)**2) :: work
 
+
   m_= m                                                                                             ! copy matrix to input (will be destroyed)
   call dsyev('N','U',size(m,1),m_,size(m,1),math_eigvalsh,work,size(work,1),ierr)
   if (ierr /= 0) math_eigvalsh = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
@@ -1126,6 +1159,7 @@ pure function math_eigvalsh33(m)
   real(pReal) :: P, Q, rho, phi
   real(pReal), parameter :: TOL=1.e-14_pReal
 
+
   I = math_invariantsSym33(m)                                                                       ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206
 
   P = I(2)-I(1)**2/3.0_pReal                                                                        ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
@@ -1157,6 +1191,7 @@ pure function math_invariantsSym33(m)
   real(pReal), dimension(3,3), intent(in) :: m
   real(pReal), dimension(3) :: math_invariantsSym33
 
+
   math_invariantsSym33(1) = math_trace33(m)
   math_invariantsSym33(2) = m(1,1)*m(2,2) + m(1,1)*m(3,3) + m(2,2)*m(3,3) &
                           -(m(1,2)**2     + m(1,3)**2     + m(2,3)**2)