diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 165cdc68f..5fdd9869c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -52,8 +52,8 @@ variables:
   MPI_Intel:      "MPI/Intel/19.1.2/IntelMPI/2019"
   MPI_GNU:        "MPI/GNU/10/OpenMPI/4.1.1"
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  PETSc_Intel:    "Libraries/PETSc/3.16.0/Intel-19.1.2-IntelMPI-2019"
-  PETSc_GNU:      "Libraries/PETSc/3.16.0/GNU-10-OpenMPI-4.1.1"
+  PETSc_Intel:    "Libraries/PETSc/3.16.1/Intel-19.1.2-IntelMPI-2019"
+  PETSc_GNU:      "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.2"
   IntelMarc:      "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
diff --git a/PRIVATE b/PRIVATE
index 62be3aa39..686e7e366 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 62be3aa392eb17cd553041996c6315e47ea3b83c
+Subproject commit 686e7e3669e21a7ed92da7c989dac55b73949b05
diff --git a/examples/config/phase/W.yaml b/examples/config/phase/W.yaml
index 8be7d38a6..740fe0992 100644
--- a/examples/config/phase/W.yaml
+++ b/examples/config/phase/W.yaml
@@ -1,4 +1,4 @@
 references:
   - https://en.wikipedia.org/wiki/Tungsten
-lattice: cF
+lattice: cI
 rho: 19300.0
diff --git a/python/damask/VERSION b/python/damask/VERSION
index bd451b119..7b46fdedc 100644
--- a/python/damask/VERSION
+++ b/python/damask/VERSION
@@ -1 +1 @@
-v3.0.0-alpha5-64-g8e08af31e
+v3.0.0-alpha5-72-g8138d3a18
diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index e16001feb..69746bbe3 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -75,7 +75,7 @@ subroutine discretization_grid_init(restart)
   print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
 
 
-  if(worldrank == 0) then
+  if (worldrank == 0) then
     fileContent = IO_read(interface_geomFile)
     call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
     fname = interface_geomFile
@@ -85,7 +85,7 @@ subroutine discretization_grid_init(restart)
     call results_closeJobFile
   else
     allocate(materialAt_global(0))                                                                  ! needed for IntelMPI
-  endif
+  end if
 
 
   call MPI_Bcast(grid,3,MPI_INTEGER,0,MPI_COMM_WORLD, ierr)
@@ -100,7 +100,7 @@ subroutine discretization_grid_init(restart)
   print'(a,3(es12.5))',   ' size   x y z: ', geomSize
   print'(a,3(es12.5))',   ' origin x y z: ', origin
 
-  if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
+  if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
 
   call fftw_mpi_init
   devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), &
@@ -109,7 +109,7 @@ subroutine discretization_grid_init(restart)
                                    PETSC_COMM_WORLD, &
                                    z, &                                                             ! domain grid size along z
                                    z_offset)                                                        ! domain grid offset along z
-  if(z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes')
+  if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes')
 
   grid3       = int(z)
   grid3Offset = int(z_offset)
@@ -136,14 +136,14 @@ subroutine discretization_grid_init(restart)
 
 !--------------------------------------------------------------------------------------------------
 ! store geometry information for post processing
-  if(.not. restart) then
+  if (.not. restart) then
     call results_openJobFile
     call results_closeGroup(results_addGroup('geometry'))
     call results_addAttribute('cells', grid,    '/geometry')
     call results_addAttribute('size',  geomSize,'/geometry')
     call results_addAttribute('origin',origin,  '/geometry')
     call results_closeJobFile
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! geometry information required by the nonlocal CP model
@@ -202,20 +202,20 @@ subroutine readVTI(grid,geomSize,origin,material, &
     if (endPos < startPos) endPos = len(fileContent,kind=pI64)                                      ! end of file without new line
 
     if (.not. inFile) then
-      if(index(fileContent(startPos:endPos),'<VTKFile',kind=pI64) /= 0_pI64) then
+      if (index(fileContent(startPos:endPos),'<VTKFile',kind=pI64) /= 0_pI64) then
         inFile = .true.
         if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format')
         headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
         compressed  = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
-      endif
+      end if
     else
       if (.not. inImage) then
         if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
           inImage = .true.
           call cellsSizeOrigin(grid,geomSize,origin,fileContent(startPos:endPos))
-        endif
+        end if
       else
-        if (index(fileContent(startPos:endPos),'<CellData>',kind=pI64) /= 0_pI64) then
+        if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
           gotCellData = .true.
           do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
             if (index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. &
@@ -230,25 +230,25 @@ subroutine readVTI(grid,geomSize,origin,material, &
               s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64   ! start (no leading whitespace)
               material = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
               exit
-            endif
+            end if
             startPos = endPos + 2_pI64
             endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
-          enddo
-        endif
-      endif
-    endif
+          end do
+        end if
+      end if
+    end if
 
     if (gotCellData) exit
     startPos = endPos + 2_pI64
 
-  enddo
+  end do
 
-  if(.not. allocated(material))       call IO_error(error_ID = 844, ext_msg='material data not found')
-  if(size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
-  if(any(geomSize<=0))                call IO_error(error_ID = 844, ext_msg='size')
-  if(any(grid<1))                     call IO_error(error_ID = 844, ext_msg='grid')
+  if (.not. allocated(material))       call IO_error(error_ID = 844, ext_msg='material data not found')
+  if (size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
+  if (any(geomSize<=0))                call IO_error(error_ID = 844, ext_msg='size')
+  if (any(grid<1))                     call IO_error(error_ID = 844, ext_msg='grid')
   material = material + 1
-  if(any(material<1))                 call IO_error(error_ID = 844, ext_msg='material ID < 0')
+  if (any(material<1))                 call IO_error(error_ID = 844, ext_msg='material ID < 0')
 
   contains
 
@@ -352,11 +352,11 @@ subroutine readVTI(grid,geomSize,origin,material, &
 
     integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes
 
-    if(compressed) then
+    if (compressed) then
       bytes = asBytes_compressed(base64_str,headerType)
     else
       bytes = asBytes_uncompressed(base64_str,headerType)
-    endif
+    end if
 
   end function asBytes
 
@@ -379,17 +379,18 @@ subroutine readVTI(grid,geomSize,origin,material, &
     integer(pI64), dimension(:), allocatable :: temp, size_inflated, size_deflated
     integer(pI64) :: headerLen, nBlock, b,s,e
 
-    if    (headerType == 'UInt32') then
+
+    if      (headerType == 'UInt32') then
       temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(4_pI64)))),pI64)
       nBlock = int(temp(1),pI64)
       headerLen = 4_pI64 * (3_pI64 + nBlock)
       temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
-    elseif(headerType == 'UInt64') then
+    else if (headerType == 'UInt64') then
       temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64)
       nBlock = int(temp(1),pI64)
       headerLen = 8_pI64 * (3_pI64 + nBlock)
       temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
-    endif
+    end if
 
     allocate(size_inflated(nBlock),source=temp(2))
     size_inflated(nBlock) = merge(temp(3),temp(2),temp(3)/=0_pI64)
@@ -402,7 +403,7 @@ subroutine readVTI(grid,geomSize,origin,material, &
       s = e + 1_pI64
       e = s + size_deflated(b) - 1_pI64
       bytes(sum(size_inflated(:b-1))+1_pI64:sum(size_inflated(:b))) = zlib_inflate(bytes_inflated(s:e),size_inflated(b))
-    enddo
+    end do
 
   end function asBytes_compressed
 
@@ -424,19 +425,19 @@ subroutine readVTI(grid,geomSize,origin,material, &
     allocate(bytes(0))
 
     s=0_pI64
-    if    (headerType == 'UInt32') then
+    if      (headerType == 'UInt32') then
       do while(s+base64_nChar(4_pI64)<(len(base64_str,pI64)))
         nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64)
         bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)]
         s = s + base64_nChar(4_pI64+nByte(1))
-      enddo
-    elseif(headerType == 'UInt64') then
+      end do
+    else if (headerType == 'UInt64') then
       do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64)))
         nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64)
         bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)]
         s = s + base64_nChar(8_pI64+nByte(1))
-      enddo
-    endif
+      end do
+    end if
 
   end function asBytes_uncompressed
 
@@ -455,11 +456,11 @@ subroutine readVTI(grid,geomSize,origin,material, &
 #endif
 
     s = index(line," "//key,back=.true.)
-    if(s==0) then
+    if (s==0) then
       getXMLValue = ''
     else
       e = s + 1 + scan(line(s+1:),"'"//'"')
-      if(scan(line(s:e-2),'=') == 0) then
+      if (scan(line(s:e-2),'=') == 0) then
         getXMLValue = ''
       else
         s = e
@@ -471,8 +472,8 @@ subroutine readVTI(grid,geomSize,origin,material, &
         e = s + index(line(s:),merge("'",'"',line(s-1:s-1)=="'")) - 1
 #endif
         getXMLValue = line(s:e-1)
-      endif
-    endif
+      end if
+    end if
 
   end function
 
@@ -510,11 +511,12 @@ function IPcoordinates0(grid,geomSize,grid3Offset)
     a,b,c, &
     i
 
+
   i = 0
   do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1)
     i = i + 1
     IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal)
-  enddo; enddo; enddo
+  end do; end do; end do
 
 end function IPcoordinates0
 
@@ -538,7 +540,7 @@ pure function nodes0(grid,geomSize,grid3Offset)
   do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1)
     n = n + 1
     nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal)
-  enddo; enddo; enddo
+  end do; end do; end do
 
 end function nodes0
 
@@ -553,6 +555,7 @@ pure function cellSurfaceArea(geomSize,grid)
 
   real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea
 
+
   cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
   cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
   cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
@@ -631,7 +634,7 @@ pure function IPneighborhood(grid)
     IPneighborhood(3,5,1,e) = 6
     IPneighborhood(3,6,1,e) = 5
 
-  enddo; enddo; enddo
+  end do; end do; end do
 
 end function IPneighborhood
 
diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90
index a3856ccaa..fde8831db 100644
--- a/src/mesh/FEM_utilities.f90
+++ b/src/mesh/FEM_utilities.f90
@@ -123,9 +123,7 @@ subroutine FEM_utilities_init
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls &
                                &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew &
                                &-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 &
-                               &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25 &
-                               &-mechanical_pc_type ml -mechanical_mg_levels_ksp_type chebyshev &
-                               &-mechanical_mg_levels_pc_type sor -mechanical_pc_ml_nullspace user',ierr)
+                               &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25', ierr)
   CHKERRQ(ierr)
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('PETSc_options',defaultVal=''),ierr)
   CHKERRQ(ierr)
@@ -151,6 +149,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
 
   PetscErrorCode :: ierr
 
+
   print'(/,a)', ' ... evaluating constitutive response ......................................'
 
   call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
@@ -178,6 +177,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
   PetscScalar          :: BCValue,BCDotValue,timeinc
   PetscErrorCode       :: ierr
 
+
   call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
   call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
   if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
@@ -189,8 +189,8 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
     CHKERRQ(ierr)
     do dof = offset+comp+1, offset+numDof, numComp
       localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
-    enddo
-  enddo
+    end do
+  end do
   call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
   call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
   call VecAssemblyEnd  (localVec, ierr); CHKERRQ(ierr)
diff --git a/src/prec.f90 b/src/prec.f90
index d6d161a94..063a0d565 100644
--- a/src/prec.f90
+++ b/src/prec.f90
@@ -273,27 +273,27 @@ subroutine selfTest
 
 
   realloc_lhs_test = [1,2]
-  if (any(realloc_lhs_test/=[1,2]))        error stop 'LHS allocation'
+  if (any(realloc_lhs_test/=[1,2]))         error stop 'LHS allocation'
 
   call random_number(r)
   r = r/minval(r)
-  if(.not. all(dEq(r,r+PREAL_EPSILON)))    error stop 'dEq'
-  if(dEq(r(1),r(2)) .and. dNeq(r(1),r(2))) error stop 'dNeq'
-  if(.not. all(dEq0(r-(r+PREAL_MIN))))     error stop 'dEq0'
+  if (.not. all(dEq(r,r+PREAL_EPSILON)))    error stop 'dEq'
+  if (dEq(r(1),r(2)) .and. dNeq(r(1),r(2))) error stop 'dNeq'
+  if (.not. all(dEq0(r-(r+PREAL_MIN))))     error stop 'dEq0'
 
   ! https://www.binaryconvert.com
   ! https://www.rapidtables.com/convert/number/binary-to-decimal.html
   f = real(prec_bytesToC_FLOAT(int([-65,+11,-102,+75],C_SIGNED_CHAR)),pReal)
-  if(dNeq(f(1),20191102.0_pReal,0.0_pReal)) error stop 'prec_bytesToC_FLOAT'
+  if (dNeq(f(1),20191102.0_pReal,0.0_pReal)) error stop 'prec_bytesToC_FLOAT'
 
   f = real(prec_bytesToC_DOUBLE(int([0,0,0,-32,+119,+65,+115,65],C_SIGNED_CHAR)),pReal)
-  if(dNeq(f(1),20191102.0_pReal,0.0_pReal)) error stop 'prec_bytesToC_DOUBLE'
+  if (dNeq(f(1),20191102.0_pReal,0.0_pReal)) error stop 'prec_bytesToC_DOUBLE'
 
   i = int(prec_bytesToC_INT32_T(int([+126,+23,+52,+1],C_SIGNED_CHAR)),pInt)
-  if(i(1) /= 20191102_pInt)                 error stop 'prec_bytesToC_INT32_T'
+  if (i(1) /= 20191102_pInt)                 error stop 'prec_bytesToC_INT32_T'
 
   i = int(prec_bytesToC_INT64_T(int([+126,+23,+52,+1,0,0,0,0],C_SIGNED_CHAR)),pInt)
-  if(i(1) /= 20191102_pInt)                 error stop 'prec_bytesToC_INT64_T'
+  if (i(1) /= 20191102_pInt)                 error stop 'prec_bytesToC_INT64_T'
 
 end subroutine selfTest