diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6da114938..f3e74e030 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -46,7 +46,7 @@ variables:
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   PETSC_GNU:       "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
   PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
-  PETSC_INTEL:    "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
+  PETSC_INTEL:    "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.3.1"
   IntelMarc:      "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
@@ -79,63 +79,63 @@ mypy:
 
 
 ###################################################################################################
-test_grid_GNU:
+grid_GNU:
   stage: compile
   script:
     - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU
 
-test_mesh_GNU:
+mesh_GNU:
   stage: compile
   script:
     - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU
 
-test_grid_GNU-64bit:
+grid_GNU-64bit:
   stage: compile
   script:
     - module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit
 
-test_mesh_GNU-64bit:
+mesh_GNU-64bit:
   stage: compile
   script:
     - module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit
 
-test_grid_IntelLLVM:
+grid_IntelLLVM:
   stage: compile
   script:
     - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_IntelLLVM
 
-test_mesh_IntelLLVM:
+mesh_IntelLLVM:
   stage: compile
   script:
     - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM
 
-test_grid_Intel:
+grid_Intel:
   stage: compile
   script:
     - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
 
-test_mesh_Intel:
+mesh_Intel:
   stage: compile
   script:
     - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
 
-test_Marc:
+Marc_Intel:
   stage: compile
   script:
     - module load $IntelMarc $HDF5Marc $MSC
@@ -158,7 +158,7 @@ setup_mesh:
     - cmake -DDAMASK_SOLVER=MESH -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
     - make -j2 all install
 
-compile_Marc:
+setup_Marc:
   stage: compile
   script:
     - module load $IntelMarc $HDF5Marc $MSC
diff --git a/PRIVATE b/PRIVATE
index b898a8b55..d570b4a1f 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53
+Subproject commit d570b4a1fac5b9b99d026d3ff9bf593615e22ce5
diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
index 0ae04928d..8a3147604 100644
--- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
+++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
@@ -6,7 +6,7 @@ references:
 
 output: [xi_sl, xi_tw]
 
-N_sl: [3, 3, 6, 0, 6]  # basal, prism, -, 1. pyr<a>, -, 2. pyr<c+a>
+N_sl: [3, 3, 6, 0, 6]  # basal, prism, 1. pyr<a>, -, 2. pyr<c+a>
 N_tw: [6, 0, 6]        # tension, -, compression
 
 xi_0_sl:   [10.e+6,  55.e+6,  60.e+6, 0.,  60.e+6]
@@ -23,9 +23,13 @@ f_sat_sl-tw: 10.0
 h_0_sl-sl: 500.0e+6
 h_0_tw-tw:  50.0e+6
 h_0_tw-sl: 150.0e+6
-h_sl-sl:  [1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
-           1.0,  1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-           1.0,  1.0,  1.0, -1.0,  1.0,  1.0, -1.0,  1.0, +1.0,  1.0]  # unused entries are indicated by -1.0
+h_sl-sl: [+1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0, -1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  1.0,  1.0,
+          -1.0, -1.0, -1.0, -1.0,  1.0,  1.0,  1.0,  1.0,  1.0]  # unused entries are indicated by -1.0
 h_tw-tw: [+1.0,  1.0, -1.0, -1.0, -1.0, -1.0,  1.0, -1.0,  1.0,  1.0,
           -1.0,  1.0]  # unused entries are indicated by -1.0
 h_tw-sl:  [1.0, -1.0,  1.0, -1.0,  1.0, -1.0,  1.0, -1.0,  1.0, -1.0,
diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
index 1db6adc6b..890f580cc 100644
--- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
+++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
@@ -8,7 +8,7 @@ references:
     https://doi.org/10.1016/j.actamat.2017.05.015
 output: [gamma_sl]
 
-N_sl: [3, 3, 0, 12]   # basal, 1. prism, -, 1. pyr<c+a>
+N_sl: [3, 3, 0, 12]   # basal, prism, -, 1. pyr<c+a>
 n_sl: 20
 a_sl: 2.0
 dot_gamma_0_sl: 0.001
@@ -20,5 +20,8 @@ xi_inf_sl: [568.e+6, 150.e+7, 0.0, 3420.e+6]
 # L. Wang et al. :
 # xi_0_sl:   [127.e+6, 96.e+6, 0.0, 240.e+6]
 
-h_sl-sl: [+1.0,  1.0,  1.0,  1.0,  1.0,  1.0, -1.0, -1.0, -1.0, -1.0,
-          -1.0, -1.0, +1.0,  1.0, -1.0,  1.0,  1.0, -1.0,  1.0,  1.0] # unused entries are indicated by -1.0
+h_sl-sl: [+1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0, -1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0] # unused entries are indicated by -1.0
diff --git a/python/damask/VERSION b/python/damask/VERSION
index 11b9d06bd..808f9b405 100644
--- a/python/damask/VERSION
+++ b/python/damask/VERSION
@@ -1 +1 @@
-v3.0.0-alpha5-389-ga000e477c
+v3.0.0-alpha5-475-g160eb1c60
diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index 584a97e87..aa3c36b07 100644
--- a/python/damask/__init__.py
+++ b/python/damask/__init__.py
@@ -21,8 +21,8 @@ from ._rotation        import Rotation         # noqa
 from ._crystal         import Crystal          # noqa
 from ._orientation     import Orientation      # noqa
 from ._table           import Table            # noqa
-from ._vtk             import VTK              # noqa
 from ._colormap        import Colormap         # noqa
+from ._vtk             import VTK              # noqa
 from ._config          import Config           # noqa
 from ._configmaterial  import ConfigMaterial   # noqa
 from ._grid            import Grid             # noqa
diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py
index 2c7f01d94..1ccc292e6 100644
--- a/python/damask/_colormap.py
+++ b/python/damask/_colormap.py
@@ -250,7 +250,7 @@ class Colormap(mpl.colors.ListedColormap):
                np.logical_or (np.isnan(field), field == gap))                                       # mask NaN (and gap if present)
 
         l,r = (field[mask].min(),field[mask].max()) if bounds is None else \
-              np.array(bounds,float)[:2]
+              (bounds[0],bounds[1])
 
         delta,avg = r-l,0.5*abs(r+l)
 
diff --git a/python/damask/_grid.py b/python/damask/_grid.py
index 5006e1b40..94f7eefd0 100644
--- a/python/damask/_grid.py
+++ b/python/damask/_grid.py
@@ -4,7 +4,7 @@ import warnings
 import multiprocessing as mp
 from functools import partial
 import typing
-from typing import Union, Optional, TextIO, List, Sequence
+from typing import Union, Optional, TextIO, List, Sequence, Literal
 from pathlib import Path
 
 import numpy as np
@@ -70,7 +70,7 @@ class Grid:
               ])
 
 
-    def __copy__(self) -> "Grid":
+    def __copy__(self) -> 'Grid':
         """Create deep copy."""
         return copy.deepcopy(self)
 
@@ -161,7 +161,7 @@ class Grid:
 
 
     @staticmethod
-    def load(fname: Union[str, Path]) -> "Grid":
+    def load(fname: Union[str, Path]) -> 'Grid':
         """
         Load from VTK image data file.
 
@@ -190,7 +190,7 @@ class Grid:
 
     @typing. no_type_check
     @staticmethod
-    def load_ASCII(fname)-> "Grid":
+    def load_ASCII(fname)-> 'Grid':
         """
         Load from geom file.
 
@@ -264,7 +264,7 @@ class Grid:
 
 
     @staticmethod
-    def load_Neper(fname: Union[str, Path]) -> "Grid":
+    def load_Neper(fname: Union[str, Path]) -> 'Grid':
         """
         Load from Neper VTK file.
 
@@ -279,7 +279,7 @@ class Grid:
             Grid-based geometry from file.
 
         """
-        v = VTK.load(fname,'vtkImageData')
+        v = VTK.load(fname,'ImageData')
         cells = np.array(v.vtk_data.GetDimensions())-1
         bbox  = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
 
@@ -292,7 +292,7 @@ class Grid:
     def load_DREAM3D(fname: Union[str, Path],
                      feature_IDs: str = None, cell_data: str = None,
                      phases: str = 'Phases', Euler_angles: str = 'EulerAngles',
-                     base_group: str = None) -> "Grid":
+                     base_group: str = None) -> 'Grid':
         """
         Load DREAM.3D (HDF5) file.
 
@@ -354,7 +354,7 @@ class Grid:
     @staticmethod
     def from_table(table: Table,
                    coordinates: str,
-                   labels: Union[str, Sequence[str]]) -> "Grid":
+                   labels: Union[str, Sequence[str]]) -> 'Grid':
         """
         Create grid from ASCII table.
 
@@ -422,13 +422,14 @@ class Grid:
             Grid-based geometry from tessellation.
 
         """
+        weights_p: FloatSequence
         if periodic:
             weights_p = np.tile(weights,27)                                                         # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
             seeds_p = np.vstack((seeds  -np.array([size[0],0.,0.]),seeds,  seeds  +np.array([size[0],0.,0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
         else:
-            weights_p = weights
+            weights_p = np.array(weights,float)
             seeds_p   = seeds
 
         coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
@@ -452,7 +453,7 @@ class Grid:
                                   size: FloatSequence,
                                   seeds: np.ndarray,
                                   material: IntSequence = None,
-                                  periodic: bool = True) -> "Grid":
+                                  periodic: bool = True) -> 'Grid':
         """
         Create grid from Voronoi tessellation.
 
@@ -538,7 +539,7 @@ class Grid:
                              surface: str,
                              threshold: float = 0.0,
                              periods: int = 1,
-                             materials: IntSequence = (0,1)) -> "Grid":
+                             materials: IntSequence = (0,1)) -> 'Grid':
         """
         Create grid from definition of triply periodic minimal surface.
 
@@ -674,7 +675,7 @@ class Grid:
 
     def show(self) -> None:
         """Show on screen."""
-        VTK.from_rectilinear_grid(self.cells,self.size,self.origin).show()
+        VTK.from_image_data(self.cells,self.size,self.origin).show()
 
 
     def add_primitive(self,
@@ -684,7 +685,7 @@ class Grid:
                       fill: int = None,
                       R: Rotation = Rotation(),
                       inverse: bool = False,
-                      periodic: bool = True) -> "Grid":
+                      periodic: bool = True) -> 'Grid':
         """
         Insert a primitive geometric object at a given position.
 
@@ -769,7 +770,7 @@ class Grid:
                    )
 
 
-    def mirror(self, directions: Sequence[str], reflect: bool = False) -> "Grid":
+    def mirror(self, directions: Sequence[str], reflect: bool = False) -> 'Grid':
         """
         Mirror grid along given directions.
 
@@ -821,7 +822,7 @@ class Grid:
                    )
 
 
-    def flip(self, directions: Sequence[str]) -> "Grid":
+    def flip(self, directions: Union[Literal['x', 'y', 'z'], Sequence[Literal['x', 'y', 'z']]]) -> 'Grid':
         """
         Flip grid along given directions.
 
@@ -841,7 +842,8 @@ class Grid:
         if not set(directions).issubset(valid):
             raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
 
-        mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
+
+        mat = np.flip(self.material, [valid.index(d) for d in directions if d in valid])
 
         return Grid(material = mat,
                     size     = self.size,
@@ -850,7 +852,7 @@ class Grid:
                    )
 
 
-    def scale(self, cells: IntSequence, periodic: bool = True) -> "Grid":
+    def scale(self, cells: IntSequence, periodic: bool = True) -> 'Grid':
         """
         Scale grid to new cells.
 
@@ -897,7 +899,7 @@ class Grid:
     def clean(self,
               stencil: int = 3,
               selection: IntSequence = None,
-              periodic: bool = True) -> "Grid":
+              periodic: bool = True) -> 'Grid':
         """
         Smooth grid by selecting most frequent material index within given stencil at each location.
 
@@ -937,7 +939,7 @@ class Grid:
                    )
 
 
-    def renumber(self) -> "Grid":
+    def renumber(self) -> 'Grid':
         """
         Renumber sorted material indices as 0,...,N-1.
 
@@ -956,7 +958,7 @@ class Grid:
                    )
 
 
-    def rotate(self, R: Rotation, fill: int = None) -> "Grid":
+    def rotate(self, R: Rotation, fill: int = None) -> 'Grid':
         """
         Rotate grid (pad if required).
 
@@ -973,14 +975,13 @@ class Grid:
             Updated grid-based geometry.
 
         """
-        if fill is None: fill = np.nanmax(self.material) + 1
-        dtype = float if isinstance(fill,float) or self.material.dtype in np.sctypes['float'] else int
-
         material = self.material
         # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
         # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
         for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]):
-            material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,output=dtype,cval=fill)
+            material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,
+                                           output=self.material.dtype,
+                                           cval=np.nanmax(self.material) + 1 if fill is None else fill)
             # avoid scipy interpolation errors for rotations close to multiples of 90°
             material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
                        np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
@@ -997,7 +998,7 @@ class Grid:
     def canvas(self,
                cells: IntSequence = None,
                offset: IntSequence = None,
-               fill: int = None) -> "Grid":
+               fill: int = None) -> 'Grid':
         """
         Crop or enlarge/pad grid.
 
@@ -1031,10 +1032,8 @@ class Grid:
         """
         offset_ = np.array(offset,int) if offset is not None else np.zeros(3,int)
         cells_ = np.array(cells,int) if cells is not None else self.cells
-        if fill is None: fill = np.nanmax(self.material) + 1
-        dtype = float if int(fill) != fill or self.material.dtype in np.sctypes['float'] else int
 
-        canvas = np.full(cells_,fill,dtype)
+        canvas = np.full(cells_,np.nanmax(self.material) + 1 if fill is None else fill,self.material.dtype)
 
         LL = np.clip( offset_,           0,np.minimum(self.cells,     cells_+offset_))
         UR = np.clip( offset_+cells_,    0,np.minimum(self.cells,     cells_+offset_))
@@ -1050,7 +1049,7 @@ class Grid:
                    )
 
 
-    def substitute(self, from_material: IntSequence, to_material: IntSequence) -> "Grid":
+    def substitute(self, from_material: IntSequence, to_material: IntSequence) -> 'Grid':
         """
         Substitute material indices.
 
@@ -1067,20 +1066,18 @@ class Grid:
             Updated grid-based geometry.
 
         """
-        def mp(entry, mapper):
-            return mapper[entry] if entry in mapper else entry
+        material = self.material.copy()
+        for f,t in zip(from_material,to_material):        # ToDo Python 3.10 has strict mode for zip
+            material[self.material==f] = t
 
-        mp = np.vectorize(mp)
-        mapper = dict(zip(from_material,to_material))
-
-        return Grid(material = mp(self.material,mapper).reshape(self.cells),
+        return Grid(material = material,
                     size     = self.size,
                     origin   = self.origin,
                     comments = self.comments+[util.execution_stamp('Grid','substitute')],
                    )
 
 
-    def sort(self) -> "Grid":
+    def sort(self) -> 'Grid':
         """
         Sort material indices such that min(material) is located at (0,0,0).
 
@@ -1106,7 +1103,7 @@ class Grid:
                         vicinity: int = 1,
                         offset: int = None,
                         trigger: IntSequence = [],
-                        periodic: bool = True) -> "Grid":
+                        periodic: bool = True) -> 'Grid':
         """
         Offset material index of points in the vicinity of xxx.
 
diff --git a/python/damask/_result.py b/python/damask/_result.py
index c87b51a89..02d4174c9 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -981,7 +981,7 @@ class Result:
             t = 'tensor'
             if o is None: o = 'fro'
         else:
-            raise ValueError(f'invalid norm order {ord}')
+            raise ValueError(f'invalid shape of {x["label"]}')
 
         return {
                 'data':  np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index fd19ec31b..8e5f7f7f9 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -812,8 +812,7 @@ class Rotation:
         return Rotation.from_basis(R)
 
     @staticmethod
-    def from_parallel(a,b,
-                      **kwargs):
+    def from_parallel(a,b):
         """
         Initialize from pairs of two orthogonal lattice basis vectors.
 
@@ -963,8 +962,7 @@ class Rotation:
                  N = 500,
                  degrees = True,
                  fractions = True,
-                 rng_seed = None,
-                 **kwargs):
+                 rng_seed = None):
         """
         Sample discrete values from a binned orientation distribution function (ODF).
 
diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index dbbfb1b10..4e3f27e0e 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -2,6 +2,7 @@ import os
 import warnings
 import multiprocessing as mp
 from pathlib import Path
+from typing import Union, Literal, List
 
 import numpy as np
 import vtk
@@ -9,6 +10,7 @@ from vtk.util.numpy_support import numpy_to_vtk            as np_to_vtk
 from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
 from vtk.util.numpy_support import vtk_to_numpy            as vtk_to_np
 
+from ._typehints import FloatSequence, IntSequence
 from . import util
 from . import Table
 
@@ -20,7 +22,7 @@ class VTK:
     High-level interface to VTK.
     """
 
-    def __init__(self,vtk_data):
+    def __init__(self, vtk_data: vtk.vtkDataSet):
         """
         New spatial visualization.
 
@@ -36,7 +38,7 @@ class VTK:
 
 
     @staticmethod
-    def from_image_data(cells,size,origin=np.zeros(3)):
+    def from_image_data(cells: IntSequence, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
         """
         Create VTK of type vtk.vtkImageData.
 
@@ -60,13 +62,13 @@ class VTK:
         vtk_data = vtk.vtkImageData()
         vtk_data.SetDimensions(*(np.array(cells)+1))
         vtk_data.SetOrigin(*(np.array(origin)))
-        vtk_data.SetSpacing(*(size/cells))
+        vtk_data.SetSpacing(*(np.array(size)/np.array(cells)))
 
         return VTK(vtk_data)
 
 
     @staticmethod
-    def from_rectilinear_grid(grid,size,origin=np.zeros(3)):
+    def from_rectilinear_grid(grid: np.ndarray, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
         """
         Create VTK of type vtk.vtkRectilinearGrid.
 
@@ -98,7 +100,7 @@ class VTK:
 
 
     @staticmethod
-    def from_unstructured_grid(nodes,connectivity,cell_type):
+    def from_unstructured_grid(nodes: np.ndarray, connectivity: np.ndarray, cell_type: str) -> 'VTK':
         """
         Create VTK of type vtk.vtkUnstructuredGrid.
 
@@ -138,7 +140,7 @@ class VTK:
 
 
     @staticmethod
-    def from_poly_data(points):
+    def from_poly_data(points: np.ndarray) -> 'VTK':
         """
         Create VTK of type vtk.polyData.
 
@@ -172,15 +174,17 @@ class VTK:
 
 
     @staticmethod
-    def load(fname,dataset_type=None):
+    def load(fname: Union[str, Path],
+             dataset_type: Literal['ImageData', 'UnstructuredGrid', 'PolyData'] = None) -> 'VTK':
         """
         Load from VTK file.
 
         Parameters
         ----------
         fname : str or pathlib.Path
-            Filename for reading. Valid extensions are .vti, .vtr, .vtu, .vtp, and .vtk.
-        dataset_type : {'vtkImageData', ''vtkRectilinearGrid', 'vtkUnstructuredGrid', 'vtkPolyData'}, optional
+            Filename for reading.
+            Valid extensions are .vti, .vtr, .vtu, .vtp, and .vtk.
+        dataset_type : {'ImageData', 'UnstructuredGrid', 'PolyData'}, optional
             Name of the vtk.vtkDataSet subclass when opening a .vtk file.
 
         Returns
@@ -234,7 +238,7 @@ class VTK:
     def _write(writer):
         """Wrapper for parallel writing."""
         writer.Write()
-    def save(self,fname,parallel=True,compress=True):
+    def save(self, fname: Union[str, Path], parallel: bool = True, compress: bool = True):
         """
         Save as VTK file.
 
@@ -280,7 +284,7 @@ class VTK:
 
     # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
     # Needs support for damask.Table
-    def add(self,data,label=None):
+    def add(self, data: Union[np.ndarray, np.ma.MaskedArray], label: str = None):
         """
         Add data to either cells or points.
 
@@ -327,7 +331,7 @@ class VTK:
             raise TypeError
 
 
-    def get(self,label):
+    def get(self, label: str) -> np.ndarray:
         """
         Get either cell or point data.
 
@@ -369,7 +373,7 @@ class VTK:
             raise ValueError(f'Array "{label}" not found.')
 
 
-    def get_comments(self):
+    def get_comments(self) -> List[str]:
         """Return the comments."""
         fielddata = self.vtk_data.GetFieldData()
         for a in range(fielddata.GetNumberOfArrays()):
@@ -379,7 +383,7 @@ class VTK:
         return []
 
 
-    def set_comments(self,comments):
+    def set_comments(self, comments: Union[str, List[str]]):
         """
         Set comments.
 
@@ -396,7 +400,7 @@ class VTK:
         self.vtk_data.GetFieldData().AddArray(s)
 
 
-    def add_comments(self,comments):
+    def add_comments(self, comments: Union[str, List[str]]):
         """
         Add comments.
 
@@ -409,7 +413,7 @@ class VTK:
         self.set_comments(self.get_comments() + ([comments] if isinstance(comments,str) else comments))
 
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """ASCII representation of the VTK data."""
         writer = vtk.vtkDataSetWriter()
         writer.SetHeader(f'# {util.execution_stamp("VTK")}')
@@ -419,7 +423,7 @@ class VTK:
         return writer.GetOutputString()
 
 
-    def show(self) -> None:
+    def show(self):
         """
         Render.
 
@@ -442,22 +446,24 @@ class VTK:
 
         mapper = vtk.vtkDataSetMapper()
         mapper.SetInputData(self.vtk_data)
+
         actor = vtk.vtkActor()
         actor.SetMapper(mapper)
+        actor.GetProperty().SetColor(230/255,150/255,68/255)
 
         ren = vtk.vtkRenderer()
+        ren.AddActor(actor)
+        ren.SetBackground(67/255,128/255,208/255)
 
         window = vtk.vtkRenderWindow()
         window.AddRenderer(ren)
-
-        ren.AddActor(actor)
-        ren.SetBackground(0.2,0.2,0.2)
-
         window.SetSize(width,height)
+        window.SetWindowName(util.execution_stamp('VTK','show'))
 
         iren = vtk.vtkRenderWindowInteractor()
         iren.SetRenderWindow(window)
-
-        iren.Initialize()
-        window.Render()
-        iren.Start()
+        if os.name == 'posix' and 'DISPLAY' not in os.environ:
+            print('Found no rendering device')
+        else:
+            window.Render()
+            iren.Start()
diff --git a/python/damask/seeds.py b/python/damask/seeds.py
index 6ff150df3..e6d1f9613 100644
--- a/python/damask/seeds.py
+++ b/python/damask/seeds.py
@@ -79,7 +79,7 @@ def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, dis
 
     s = 1
     i = 0
-    progress = _util._ProgressBar(N_seeds+1,'',50)
+    progress = _util.ProgressBar(N_seeds+1,'',50)
     while s < N_seeds:
         i += 1
         candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
diff --git a/python/damask/util.py b/python/damask/util.py
index 0581302db..2872762b9 100644
--- a/python/damask/util.py
+++ b/python/damask/util.py
@@ -7,12 +7,16 @@ import subprocess
 import shlex
 import re
 import fractions
+import collections.abc as abc
 from functools import reduce
+from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, Optional
+from pathlib import Path
 
 import numpy as np
 import h5py
 
 from . import version
+from ._typehints import FloatSequence
 
 # limit visibility
 __all__=[
@@ -50,16 +54,16 @@ _colors = {
 ####################################################################################################
 # Functions
 ####################################################################################################
-def srepr(arg,glue = '\n'):
+def srepr(msg, glue: str = '\n') -> str:
     r"""
     Join items with glue string.
 
     Parameters
     ----------
-    arg : iterable
+    msg : object with __repr__ or sequence of objects with __repr__
         Items to join.
     glue : str, optional
-        Glue used for joining operation. Defaults to \n.
+        Glue used for joining operation. Defaults to '\n'.
 
     Returns
     -------
@@ -67,21 +71,21 @@ def srepr(arg,glue = '\n'):
         String representation of the joined items.
 
     """
-    if (not hasattr(arg, 'strip') and
-           (hasattr(arg, '__getitem__') or
-            hasattr(arg, '__iter__'))):
-        return glue.join(str(x) for x in arg)
+    if (not hasattr(msg, 'strip') and
+           (hasattr(msg, '__getitem__') or
+            hasattr(msg, '__iter__'))):
+        return glue.join(str(x) for x in msg)
     else:
-       return arg if isinstance(arg,str) else repr(arg)
+       return msg if isinstance(msg,str) else repr(msg)
 
 
-def emph(what):
+def emph(msg) -> str:
     """
     Format with emphasis.
 
     Parameters
     ----------
-    what : object with __repr__ or iterable of objects with __repr__.
+    msg : object with __repr__ or sequence of objects with __repr__
         Message to format.
 
     Returns
@@ -90,15 +94,15 @@ def emph(what):
         Formatted string representation of the joined items.
 
     """
-    return _colors['bold']+srepr(what)+_colors['end_color']
+    return _colors['bold']+srepr(msg)+_colors['end_color']
 
-def deemph(what):
+def deemph(msg) -> str:
     """
     Format with deemphasis.
 
     Parameters
     ----------
-    what : object with __repr__ or iterable of objects with __repr__.
+    msg : object with __repr__ or sequence of objects with __repr__
         Message to format.
 
     Returns
@@ -107,15 +111,15 @@ def deemph(what):
         Formatted string representation of the joined items.
 
     """
-    return _colors['dim']+srepr(what)+_colors['end_color']
+    return _colors['dim']+srepr(msg)+_colors['end_color']
 
-def warn(what):
+def warn(msg) -> str:
     """
     Format for warning.
 
     Parameters
     ----------
-    what : object with __repr__ or iterable of objects with __repr__.
+    msg : object with __repr__ or sequence of objects with __repr__
         Message to format.
 
     Returns
@@ -124,15 +128,15 @@ def warn(what):
         Formatted string representation of the joined items.
 
     """
-    return _colors['warning']+emph(what)+_colors['end_color']
+    return _colors['warning']+emph(msg)+_colors['end_color']
 
-def strikeout(what):
+def strikeout(msg) -> str:
     """
     Format as strikeout.
 
     Parameters
     ----------
-    what : object with __repr__ or iterable of objects with __repr__.
+    msg : object with __repr__ or iterable of objects with __repr__
         Message to format.
 
     Returns
@@ -141,10 +145,10 @@ def strikeout(what):
         Formatted string representation of the joined items.
 
     """
-    return _colors['crossout']+srepr(what)+_colors['end_color']
+    return _colors['crossout']+srepr(msg)+_colors['end_color']
 
 
-def run(cmd,wd='./',env=None,timeout=None):
+def run(cmd: str, wd: str = './', env: Dict[str, str] = None, timeout: int = None) -> Tuple[str, str]:
     """
     Run a command.
 
@@ -153,7 +157,7 @@ def run(cmd,wd='./',env=None,timeout=None):
     cmd : str
         Command to be executed.
     wd : str, optional
-        Working directory of process. Defaults to ./ .
+        Working directory of process. Defaults to './'.
     env : dict, optional
         Environment for execution.
     timeout : integer, optional
@@ -185,7 +189,7 @@ def run(cmd,wd='./',env=None,timeout=None):
 execute = run
 
 
-def natural_sort(key):
+def natural_sort(key: str) -> List[Union[int, str]]:
     """
     Natural sort.
 
@@ -200,7 +204,10 @@ def natural_sort(key):
     return [ convert(c) for c in re.split('([0-9]+)', key) ]
 
 
-def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
+def show_progress(iterable: Iterable,
+                  N_iter: int = None,
+                  prefix: str = '',
+                  bar_length: int = 50) -> Any:
     """
     Decorate a loop with a progress bar.
 
@@ -208,39 +215,49 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
 
     Parameters
     ----------
-    iterable : iterable or function with yield statement
-        Iterable (or function with yield statement) to be decorated.
+    iterable : iterable
+        Iterable to be decorated.
     N_iter : int, optional
-        Total number of iterations. Required unless obtainable as len(iterable).
+        Total number of iterations. Required if iterable is not a sequence.
     prefix : str, optional
         Prefix string.
     bar_length : int, optional
         Length of progress bar in characters. Defaults to 50.
 
     """
-    if N_iter in [0,1] or (hasattr(iterable,'__len__') and len(iterable) <= 1):
+    if isinstance(iterable,abc.Sequence):
+       if N_iter is None:
+          N = len(iterable)
+       else:
+          raise ValueError('N_iter given for sequence')
+    else:
+       if N_iter is None:
+          raise ValueError('N_iter not given')
+       else:
+          N = N_iter
+
+    if N <= 1:
         for item in iterable:
             yield item
     else:
-        status = _ProgressBar(N_iter if N_iter is not None else len(iterable),prefix,bar_length)
-
+        status = ProgressBar(N,prefix,bar_length)
         for i,item in enumerate(iterable):
             yield item
             status.update(i)
 
 
-def scale_to_coprime(v):
+def scale_to_coprime(v: FloatSequence) -> np.ndarray:
     """
     Scale vector to co-prime (relatively prime) integers.
 
     Parameters
     ----------
-    v : numpy.ndarray of shape (:)
+    v : sequence of float, len (:)
         Vector to scale.
 
     Returns
     -------
-    m : numpy.ndarray of shape (:)
+    m : numpy.ndarray, shape (:)
         Vector scaled to co-prime numbers.
 
     """
@@ -257,17 +274,21 @@ def scale_to_coprime(v):
         except AttributeError:
             return a * b // np.gcd(a, b)
 
-    m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(int)
+    v_ = np.array(v)
+    m = (v_ * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v_))**0.5).astype(int)
     m = m//reduce(np.gcd,m)
 
     with np.errstate(invalid='ignore'):
-        if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]):
-            raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?')
+        if not np.allclose(np.ma.masked_invalid(v_/m),v_[np.argmax(abs(v_))]/m[np.argmax(abs(v_))]):
+            raise ValueError(f'Invalid result {m} for input {v_}. Insufficient precision?')
 
     return m
 
 
-def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
+def project_equal_angle(vector: np.ndarray,
+                        direction: Literal['x', 'y', 'z'] = 'z',
+                        normalize: bool = True,
+                        keepdims: bool = False) -> np.ndarray:
     """
     Apply equal-angle projection to vector.
 
@@ -275,9 +296,8 @@ def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
     ----------
     vector : numpy.ndarray, shape (...,3)
         Vector coordinates to be projected.
-    direction : str
-        Projection direction 'x', 'y', or 'z'.
-        Defaults to 'z'.
+    direction : {'x', 'y', 'z'}
+        Projection direction. Defaults to 'z'.
     normalize : bool
         Ensure unit length of input vector. Defaults to True.
     keepdims : bool
@@ -309,7 +329,10 @@ def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
     return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
                    -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
-def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
+def project_equal_area(vector: np.ndarray,
+                       direction: Literal['x', 'y', 'z'] = 'z',
+                       normalize: bool = True,
+                       keepdims: bool = False) -> np.ndarray:
     """
     Apply equal-area projection to vector.
 
@@ -317,9 +340,8 @@ def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
     ----------
     vector : numpy.ndarray, shape (...,3)
         Vector coordinates to be projected.
-    direction : str
-        Projection direction 'x', 'y', or 'z'.
-        Defaults to 'z'.
+    direction : {'x', 'y', 'z'}
+        Projection direction. Defaults to 'z'.
     normalize : bool
         Ensure unit length of input vector. Defaults to True.
     keepdims : bool
@@ -351,15 +373,14 @@ def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
     return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
                    -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
-
-def execution_stamp(class_name,function_name=None):
+def execution_stamp(class_name: str, function_name: str = None) -> str:
     """Timestamp the execution of a (function within a) class."""
     now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
     _function_name = '' if function_name is None else f'.{function_name}'
     return f'damask.{class_name}{_function_name} v{version} ({now})'
 
 
-def hybrid_IA(dist,N,rng_seed=None):
+def hybrid_IA(dist: np.ndarray, N: int, rng_seed = None) -> np.ndarray:
     """
     Hybrid integer approximation.
 
@@ -387,7 +408,10 @@ def hybrid_IA(dist,N,rng_seed=None):
     return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]]
 
 
-def shapeshifter(fro,to,mode='left',keep_ones=False):
+def shapeshifter(fro: Tuple[int, ...],
+                 to: Tuple[int, ...],
+                 mode: Literal['left','right'] = 'left',
+                 keep_ones: bool = False) -> Tuple[Optional[int], ...]:
     """
     Return dimensions that reshape 'fro' to become broadcastable to 'to'.
 
@@ -398,9 +422,9 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
     to : tuple
         Target shape of array after broadcasting.
         len(to) cannot be less than len(fro).
-    mode : str, optional
+    mode : {'left', 'right'}, optional
         Indicates whether new axes are preferably added to
-        either 'left' or 'right' of the original shape.
+        either left or right of the original shape.
         Defaults to 'left'.
     keep_ones : bool, optional
         Treat '1' in fro as literal value instead of dimensional placeholder.
@@ -434,21 +458,22 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
     fro = (1,) if not len(fro) else fro
     to  = (1,) if not len(to)  else to
     try:
-        grp = re.match(beg[mode]
+        match = re.match(beg[mode]
                       +f',{sep[mode]}'.join(map(lambda x: f'{x}'
                                                           if x>1 or (keep_ones and len(fro)>1) else
                                                           '\\d+',fro))
-                      +f',{end[mode]}',
-                       ','.join(map(str,to))+',').groups()
-    except AttributeError:
+                      +f',{end[mode]}',','.join(map(str,to))+',')
+        assert match
+        grp = match.groups()
+    except AssertionError:
         raise ValueError(f'Shapes can not be shifted {fro} --> {to}')
-    fill = ()
+    fill: Tuple[Optional[int], ...] = ()
     for g,d in zip(grp,fro+(None,)):
         fill += (1,)*g.count(',')+(d,)
     return fill[:-1]
 
 
-def shapeblender(a,b):
+def shapeblender(a: Tuple[int, ...], b: Tuple[int, ...]) -> Tuple[int, ...]:
     """
     Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
 
@@ -476,7 +501,7 @@ def shapeblender(a,b):
     return a + b[i:]
 
 
-def extend_docstring(extra_docstring):
+def extend_docstring(extra_docstring: str) -> Callable:
     """
     Decorator: Append to function's docstring.
 
@@ -492,7 +517,7 @@ def extend_docstring(extra_docstring):
     return _decorator
 
 
-def extended_docstring(f,extra_docstring):
+def extended_docstring(f: Callable, extra_docstring: str) -> Callable:
     """
     Decorator: Combine another function's docstring with a given docstring.
 
@@ -510,7 +535,7 @@ def extended_docstring(f,extra_docstring):
     return _decorator
 
 
-def DREAM3D_base_group(fname):
+def DREAM3D_base_group(fname: Union[str, Path]) -> str:
     """
     Determine the base group of a DREAM.3D file.
 
@@ -536,7 +561,7 @@ def DREAM3D_base_group(fname):
 
     return base_group
 
-def DREAM3D_cell_data_group(fname):
+def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
     """
     Determine the cell data group of a DREAM.3D file.
 
@@ -568,18 +593,18 @@ def DREAM3D_cell_data_group(fname):
     return cell_data_group
 
 
-def Bravais_to_Miller(*,uvtw=None,hkil=None):
+def Bravais_to_Miller(*, uvtw: np.ndarray = None, hkil: np.ndarray = None) -> np.ndarray:
     """
     Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
 
     Parameters
     ----------
-    uvtw|hkil : numpy.ndarray of shape (...,4)
+    uvtw|hkil : numpy.ndarray, shape (...,4)
         Miller–Bravais indices of crystallographic direction [uvtw] or plane normal (hkil).
 
     Returns
     -------
-    uvw|hkl : numpy.ndarray of shape (...,3)
+    uvw|hkl : numpy.ndarray, shape (...,3)
         Miller indices of [uvw] direction or (hkl) plane normal.
 
     """
@@ -595,18 +620,18 @@ def Bravais_to_Miller(*,uvtw=None,hkil=None):
     return np.einsum('il,...l',basis,axis)
 
 
-def Miller_to_Bravais(*,uvw=None,hkl=None):
+def Miller_to_Bravais(*, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray:
     """
     Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
 
     Parameters
     ----------
-    uvw|hkl : numpy.ndarray of shape (...,3)
+    uvw|hkl : numpy.ndarray, shape (...,3)
         Miller indices of crystallographic direction [uvw] or plane normal (hkl).
 
     Returns
     -------
-    uvtw|hkil : numpy.ndarray of shape (...,4)
+    uvtw|hkil : numpy.ndarray, shape (...,4)
         Miller–Bravais indices of [uvtw] direction or (hkil) plane normal.
 
     """
@@ -624,7 +649,7 @@ def Miller_to_Bravais(*,uvw=None,hkl=None):
     return np.einsum('il,...l',basis,axis)
 
 
-def dict_prune(d):
+def dict_prune(d: Dict) -> Dict:
     """
     Recursively remove empty dictionaries.
 
@@ -650,7 +675,7 @@ def dict_prune(d):
     return new
 
 
-def dict_flatten(d):
+def dict_flatten(d: Dict) -> Dict:
     """
     Recursively remove keys of single-entry dictionaries.
 
@@ -678,14 +703,14 @@ def dict_flatten(d):
 ####################################################################################################
 # Classes
 ####################################################################################################
-class _ProgressBar:
+class ProgressBar:
     """
     Report progress of an interation as a status bar.
 
     Works for 0-based loops, ETA is estimated by linear extrapolation.
     """
 
-    def __init__(self,total,prefix,bar_length):
+    def __init__(self, total: int, prefix: str, bar_length: int):
         """
         Set current time as basis for ETA estimation.
 
@@ -708,7 +733,7 @@ class _ProgressBar:
         sys.stderr.write(f"{self.prefix} {'░'*self.bar_length}   0% ETA n/a")
         sys.stderr.flush()
 
-    def update(self,iteration):
+    def update(self, iteration: int) -> None:
 
         fraction = (iteration+1) / self.total
         filled_length = int(self.bar_length * fraction)
diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py
index 156907670..2321745aa 100644
--- a/python/tests/test_Colormap.py
+++ b/python/tests/test_Colormap.py
@@ -81,6 +81,7 @@ class TestColormap:
         assert Colormap.from_predefined('strain') == Colormap.from_predefined('strain')
         assert Colormap.from_predefined('strain') != Colormap.from_predefined('stress')
         assert Colormap.from_predefined('strain',N=128) != Colormap.from_predefined('strain',N=64)
+        assert not Colormap.from_predefined('strain',N=128) == 1
 
     @pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)),
                                          ([0,0,0],[1,1,1])])
diff --git a/python/tests/test_Crystal.py b/python/tests/test_Crystal.py
index 9c6de965d..235ed554a 100644
--- a/python/tests/test_Crystal.py
+++ b/python/tests/test_Crystal.py
@@ -40,6 +40,9 @@ class TestCrystal:
                     alpha=alpha,beta=beta,gamma=gamma)
         assert np.allclose(np.eye(3),np.einsum('ik,jk',c.basis_real,c.basis_reciprocal))
 
+    def test_basis_invalid(self):
+        with pytest.raises(KeyError):
+            Crystal(family='cubic').basis_real
 
     @pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
     @pytest.mark.parametrize('vector',np.array([
diff --git a/python/tests/test_Grid.py b/python/tests/test_Grid.py
index 2b9ca22f4..6dd94f4bb 100644
--- a/python/tests/test_Grid.py
+++ b/python/tests/test_Grid.py
@@ -44,6 +44,7 @@ class TestGrid:
 
     def test_equal(self,default):
         assert default == default
+        assert not default == 42
 
     def test_repr(self,default):
         print(default)
diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py
index 9f623fc6c..2a32adea4 100644
--- a/python/tests/test_Orientation.py
+++ b/python/tests/test_Orientation.py
@@ -364,6 +364,11 @@ class TestOrientation:
                 table.save(reference)
             assert np.allclose(P,Table.load(reference).get('Schmid'))
 
+    def test_Schmid_invalid(self):
+        with pytest.raises(KeyError):
+            Orientation(lattice='fcc').Schmid()
+
+
 ### vectorization tests ###
 
     @pytest.mark.parametrize('lattice',['hP','cI','cF']) # tI not included yet
@@ -505,3 +510,7 @@ class TestOrientation:
         for loc in np.random.randint(0,blend,(10,len(blend))):
             assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]),
                                o.to_pole(uvw=v)[tuple(loc)])
+
+    def test_mul_invalid(self):
+        with pytest.raises(TypeError):
+            Orientation.from_random(lattice='cF')*np.ones(3)
diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py
index 4fea08eb3..1da56dedf 100644
--- a/python/tests/test_Result.py
+++ b/python/tests/test_Result.py
@@ -102,6 +102,9 @@ class TestResult:
         with pytest.raises(AttributeError):
             default.view('invalid',True)
 
+    def test_add_invalid(self,default):
+        default.add_absolute('xxxx')
+
     def test_add_absolute(self,default):
         default.add_absolute('F_e')
         in_memory = np.abs(default.place('F_e'))
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 2d623adf5..a431bc64b 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -792,6 +792,11 @@ class TestRotation:
         R = Rotation.from_random(shape,rng_seed=1)
         assert R == R if shape is None else (R == R).all()
 
+    @pytest.mark.parametrize('shape',[None,5,(4,6)])
+    def test_allclose(self,shape):
+        R = Rotation.from_random(shape,rng_seed=1)
+        assert R.allclose(R)
+
     @pytest.mark.parametrize('shape',[None,5,(4,6)])
     def test_unequal(self,shape):
         R = Rotation.from_random(shape,rng_seed=1)
@@ -1124,3 +1129,7 @@ class TestRotation:
         weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
 
         assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5
+
+    def test_mul_invalid(self):
+        with pytest.raises(TypeError):
+            Rotation.from_random()*np.ones(3)
diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py
index 68224ff33..e59409a20 100644
--- a/python/tests/test_VTK.py
+++ b/python/tests/test_VTK.py
@@ -28,6 +28,10 @@ class TestVTK:
     def _patch_execution_stamp(self, patch_execution_stamp):
         print('patched damask.util.execution_stamp')
 
+    def test_show(sef,default,monkeypatch):
+        monkeypatch.delenv('DISPLAY',raising=False)
+        default.show()
+
     def test_rectilinearGrid(self,tmp_path):
         cells  = np.random.randint(5,10,3)*2
         size   = np.random.random(3) + 1.0
diff --git a/src/YAML_types.f90 b/src/YAML_types.f90
index 5884c2850..c57f602cf 100644
--- a/src/YAML_types.f90
+++ b/src/YAML_types.f90
@@ -62,6 +62,8 @@ module YAML_types
       tNode_get_byKey_as1dString   => tNode_get_byKey_as1dString
     procedure :: &
       getKey                       => tNode_get_byIndex_asKey
+    procedure :: &
+      Keys                         => tNode_getKeys
     procedure :: &
       getIndex                     => tNode_get_byKey_asIndex
     procedure :: &
@@ -117,7 +119,7 @@ module YAML_types
 
   type, extends(tNode), public :: tList
 
-    class(tItem), pointer  :: first => null()
+    class(tItem), pointer  :: first => NULL()
 
     contains
     procedure :: asFormattedString => tList_asFormattedString
@@ -144,8 +146,8 @@ module YAML_types
 
   type :: tItem
     character(len=:), allocatable :: key
-    class(tNode),     pointer     :: node => null()
-    class(tItem),     pointer     :: next => null()
+    class(tNode),     pointer     :: node => NULL()
+    class(tItem),     pointer     :: next => NULL()
 
     contains
     final :: tItem_finalize
@@ -221,22 +223,22 @@ subroutine selfTest
     select type(s1)
       class is(tScalar)
         s1 = '2'
-    endselect
+    end select
 
     select type(s2)
       class is(tScalar)
         s2 = '3'
-    endselect
+    end select
 
     select type(s3)
       class is(tScalar)
         s3 = '4'
-    endselect
+    end select
 
     select type(s4)
       class is(tScalar)
         s4 = '5'
-    endselect
+    end select
 
 
     allocate(tList::l1)
@@ -249,14 +251,14 @@ subroutine selfTest
         if (any(dNeq(l1%as1dFloat(),[2.0_pReal,3.0_pReal]))) error stop 'tList_as1dFloat'
         if (n%get_asInt(1) /= 2)                             error stop 'byIndex_asInt'
         if (dNeq(n%get_asFloat(2),3.0_pReal))                error stop 'byIndex_asFloat'
-    endselect
+    end select
 
     allocate(tList::l3)
     select type(l3)
       class is(tList)
         call l3%append(s3)
         call l3%append(s4)
-    endselect
+    end select
 
     allocate(tList::l2)
     select type(l2)
@@ -332,9 +334,12 @@ function tNode_asScalar(self) result(scalar)
   class(tNode),   intent(in), target :: self
   class(tScalar), pointer            :: scalar
 
+
   select type(self)
     class is(tScalar)
       scalar => self
+    class default
+      nullify(scalar)
   end select
 
 end function tNode_asScalar
@@ -348,9 +353,12 @@ function tNode_asList(self) result(list)
   class(tNode), intent(in), target :: self
   class(tList), pointer            :: list
 
+
   select type(self)
     class is(tList)
       list => self
+    class default
+      nullify(list)
   end select
 
 end function tNode_asList
@@ -364,9 +372,12 @@ function tNode_asDict(self) result(dict)
   class(tNode), intent(in), target :: self
   class(tDict), pointer            :: dict
 
+
   select type(self)
     class is(tDict)
       dict => self
+    class default
+      nullify(dict)
   end select
 
 end function tNode_asDict
@@ -385,12 +396,13 @@ function tNode_get_byIndex(self,i) result(node)
   class(tItem),  pointer :: item
   integer :: j
 
+
   select type(self)
     class is(tList)
       self_ => self%asList()
     class default
       call IO_error(706,ext_msg='Expected list')
-  endselect
+  end select
 
   item => self_%first
 
@@ -409,15 +421,14 @@ end function tNode_get_byIndex
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_asFloat(self,i) result(nodeAsFloat)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   real(pReal) :: nodeAsFloat
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  node   => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tScalar)
       scalar => node%asScalar()
       nodeAsFloat = scalar%asFloat()
@@ -433,15 +444,15 @@ end function tNode_get_byIndex_asFloat
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_asInt(self,i) result(nodeAsInt)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   integer :: nodeAsInt
 
   class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  node   => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tScalar)
       scalar => node%asScalar()
       nodeAsInt = scalar%asInt()
@@ -457,21 +468,20 @@ end function tNode_get_byIndex_asInt
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_asBool(self,i) result(nodeAsBool)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   logical :: nodeAsBool
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  node   => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tScalar)
       scalar => node%asScalar()
       nodeAsBool = scalar%asBool()
     class default
       call IO_error(706,ext_msg='Expected scalar Boolean')
-  endselect
+  end select
 
 end function tNode_get_byIndex_asBool
 
@@ -481,21 +491,20 @@ end function tNode_get_byIndex_asBool
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_asString(self,i) result(nodeAsString)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   character(len=:), allocatable :: nodeAsString
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  node   => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tScalar)
       scalar => node%asScalar()
       nodeAsString = scalar%asString()
     class default
       call IO_error(706,ext_msg='Expected scalar string')
-  endselect
+  end select
 
 end function tNode_get_byIndex_asString
 
@@ -505,21 +514,20 @@ end function tNode_get_byIndex_asString
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_as1dFloat(self,i) result(nodeAs1dFloat)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   real(pReal), dimension(:), allocatable :: nodeAs1dFloat
 
-  class(tNode), pointer  :: node
   class(tList),  pointer :: list
 
-  node => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tList)
       list => node%asList()
       nodeAs1dFloat = list%as1dFloat()
     class default
       call IO_error(706,ext_msg='Expected list of floats')
-  endselect
+  end select
 
 end function tNode_get_byIndex_as1dFloat
 
@@ -529,21 +537,20 @@ end function tNode_get_byIndex_as1dFloat
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_as1dInt(self,i) result(nodeAs1dInt)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   integer, dimension(:), allocatable :: nodeAs1dInt
 
-  class(tNode), pointer :: node
   class(tList), pointer :: list
 
-  node => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tList)
       list => node%asList()
       nodeAs1dInt = list%as1dInt()
     class default
       call IO_error(706,ext_msg='Expected list of integers')
-  endselect
+  end select
 
 end function tNode_get_byIndex_as1dInt
 
@@ -553,21 +560,20 @@ end function tNode_get_byIndex_as1dInt
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_as1dBool(self,i) result(nodeAs1dBool)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   logical, dimension(:), allocatable :: nodeAs1dBool
 
-  class(tNode), pointer :: node
   class(tList), pointer :: list
 
-  node => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tList)
       list => node%asList()
       nodeAs1dBool = list%as1dBool()
     class default
       call IO_error(706,ext_msg='Expected list of Booleans')
-  endselect
+  end select
 
 end function tNode_get_byIndex_as1dBool
 
@@ -577,21 +583,20 @@ end function tNode_get_byIndex_as1dBool
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byIndex_as1dString(self,i) result(nodeAs1dString)
 
-  class(tNode), intent(in), target :: self
-  integer,      intent(in)         :: i
+  class(tNode), intent(in) :: self
+  integer,      intent(in) :: i
   character(len=:), allocatable, dimension(:) :: nodeAs1dString
 
-  class(tNode), pointer :: node
   type(tList),  pointer :: list
 
-  node => self%get(i)
-  select type(node)
+
+  select type(node => self%get(i))
     class is(tList)
       list => node%asList()
       nodeAs1dString = list%as1dString()
     class default
       call IO_error(706,ext_msg='Expected list of strings')
-  endselect
+  end select
 
 end function tNode_get_byIndex_as1dString
 
@@ -609,22 +614,50 @@ function tNode_get_byIndex_asKey(self,i)  result(key)
   type(tDict), pointer :: dict
   type(tItem), pointer :: item
 
+
   select type(self)
     class is(tDict)
       dict => self%asDict()
       item => dict%first
       do j = 1, min(i,dict%length)-1
         item => item%next
-      enddo
+      end do
     class default
       call IO_error(706,ext_msg='Expected dict')
-  endselect
+  end select
 
   key = item%key
 
 end function tNode_get_byIndex_asKey
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Get all keys from a dictionary
+!--------------------------------------------------------------------------------------------------
+function tNode_getKeys(self) result(keys)
+
+  class(tNode), intent(in) :: self
+  character(len=:), dimension(:), allocatable :: keys
+
+  character(len=pStringLen), dimension(:), allocatable :: temp
+  integer :: j, l
+
+
+  allocate(temp(self%length))
+  l = 0
+  do j = 1, self%length
+    temp(j) = self%getKey(j)
+    l = max(len_trim(temp(j)),l)
+  end do
+
+  allocate(character(l)::keys(self%length))
+  do j = 1, self%length
+    keys(j) = trim(temp(j))
+  end do
+
+end function tNode_getKeys
+
+
 !-------------------------------------------------------------------------------------------------
 !> @brief Checks if a given key/item is present in the dict/list
 !-------------------------------------------------------------------------------------------------
@@ -658,7 +691,7 @@ function tNode_contains(self,k)  result(exists)
       enddo
     class default
       call IO_error(706,ext_msg='Expected list or dict')
-  endselect
+  end select
 
 end function tNode_contains
 
@@ -686,7 +719,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
       self_ => self%asDict()
     class default
       call IO_error(706,ext_msg='Expected dict for key '//k)
-  endselect
+  end select
 
   j = 1
   item => self_%first
@@ -713,23 +746,22 @@ end function tNode_get_byKey
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat)
 
-  class(tNode),     intent(in), target  :: self
-  character(len=*), intent(in)          :: k
-  real(pReal),      intent(in),optional :: defaultVal
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
+  real(pReal),      intent(in), optional :: defaultVal
   real(pReal) :: nodeAsFloat
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tScalar)
         scalar => node%asScalar()
         nodeAsFloat = scalar%asFloat()
       class default
         call IO_error(706,ext_msg='Expected scalar float for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAsFloat = defaultVal
   else
@@ -744,23 +776,22 @@ end function tNode_get_byKey_asFloat
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt)
 
-  class(tNode),     intent(in), target  :: self
-  character(len=*), intent(in)          :: k
-  integer,          intent(in),optional :: defaultVal
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
+  integer,          intent(in), optional :: defaultVal
   integer :: nodeAsInt
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tScalar)
         scalar => node%asScalar()
         nodeAsInt = scalar%asInt()
       class default
         call IO_error(706,ext_msg='Expected scalar integer for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAsInt = defaultVal
   else
@@ -775,23 +806,22 @@ end function tNode_get_byKey_asInt
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool)
 
-  class(tNode),     intent(in), target  :: self
-  character(len=*), intent(in)          :: k
-  logical,          intent(in),optional :: defaultVal
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
+  logical,          intent(in), optional :: defaultVal
   logical :: nodeAsBool
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tScalar)
         scalar => node%asScalar()
         nodeAsBool = scalar%asBool()
       class default
         call IO_error(706,ext_msg='Expected scalar Boolean for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAsBool = defaultVal
   else
@@ -806,23 +836,22 @@ end function tNode_get_byKey_asBool
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString)
 
-  class(tNode),     intent(in), target  :: self
-  character(len=*), intent(in)          :: k
-  character(len=*), intent(in),optional :: defaultVal
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
+  character(len=*), intent(in), optional :: defaultVal
   character(len=:), allocatable :: nodeAsString
 
-  class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tScalar)
         scalar => node%asScalar()
         nodeAsString = scalar%asString()
       class default
         call IO_error(706,ext_msg='Expected scalar string for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAsString = defaultVal
   else
@@ -837,25 +866,24 @@ end function tNode_get_byKey_asString
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs1dFloat)
 
-  class(tNode),     intent(in), target                 :: self
-  character(len=*), intent(in)                         :: k
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
   real(pReal),      intent(in), dimension(:), optional :: defaultVal
   integer,          intent(in),               optional :: requiredSize
 
   real(pReal), dimension(:), allocatable :: nodeAs1dFloat
 
-  class(tNode), pointer :: node
   type(tList),  pointer :: list
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tList)
         list => node%asList()
         nodeAs1dFloat = list%as1dFloat()
       class default
         call IO_error(706,ext_msg='Expected 1D float array for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAs1dFloat = defaultVal
   else
@@ -874,25 +902,24 @@ end function tNode_get_byKey_as1dFloat
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat)
 
-  class(tNode),     intent(in), target                   :: self
-  character(len=*), intent(in)                           :: k
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
   real(pReal),      intent(in), dimension(:,:), optional :: defaultVal
   integer,          intent(in), dimension(2),   optional :: requiredShape
 
   real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat
 
-  class(tNode), pointer :: node
   type(tList),  pointer :: rows
 
+
   if(self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tList)
         rows => node%asList()
         nodeAs2dFloat = rows%as2dFloat()
       class default
         call IO_error(706,ext_msg='Expected 2D float array for key '//k)
-    endselect
+    end select
   elseif(present(defaultVal)) then
     nodeAs2dFloat = defaultVal
   else
@@ -911,24 +938,22 @@ end function tNode_get_byKey_as2dFloat
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_as1dInt(self,k,defaultVal,requiredSize) result(nodeAs1dInt)
 
-  class(tNode),          intent(in), target   :: self
-  character(len=*),      intent(in)           :: k
+  class(tNode),          intent(in) :: self
+  character(len=*),      intent(in) :: k
   integer, dimension(:), intent(in), optional :: defaultVal
   integer,               intent(in), optional :: requiredSize
   integer, dimension(:), allocatable :: nodeAs1dInt
 
-  class(tNode), pointer :: node
   type(tList),  pointer :: list
 
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tList)
         list => node%asList()
         nodeAs1dInt = list%as1dInt()
       class default
         call IO_error(706,ext_msg='Expected 1D integer array for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAs1dInt = defaultVal
   else
@@ -947,23 +972,22 @@ end function tNode_get_byKey_as1dInt
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_as1dBool(self,k,defaultVal) result(nodeAs1dBool)
 
-  class(tNode),          intent(in), target   :: self
-  character(len=*),      intent(in)           :: k
+  class(tNode),          intent(in) :: self
+  character(len=*),      intent(in) :: k
   logical, dimension(:), intent(in), optional :: defaultVal
   logical, dimension(:), allocatable          :: nodeAs1dBool
 
-  class(tNode), pointer :: node
   type(tList),  pointer :: list
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tList)
         list => node%asList()
         nodeAs1dBool = list%as1dBool()
       class default
         call IO_error(706,ext_msg='Expected 1D Boolean array for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAs1dBool = defaultVal
   else
@@ -978,23 +1002,22 @@ end function tNode_get_byKey_as1dBool
 !--------------------------------------------------------------------------------------------------
 function tNode_get_byKey_as1dString(self,k,defaultVal) result(nodeAs1dString)
 
-  class(tNode),     intent(in), target                 :: self
-  character(len=*), intent(in)                         :: k
+  class(tNode),     intent(in) :: self
+  character(len=*), intent(in) :: k
   character(len=*), intent(in), dimension(:), optional :: defaultVal
   character(len=:), allocatable, dimension(:)          :: nodeAs1dString
 
-  class(tNode), pointer :: node
   type(tList),  pointer :: list
 
+
   if (self%contains(k)) then
-    node => self%get(k)
-    select type(node)
+    select type(node => self%get(k))
       class is(tList)
         list => node%asList()
         nodeAs1dString = list%as1dString()
       class default
         call IO_error(706,ext_msg='Expected 1D string array for key '//k)
-    endselect
+    end select
   elseif (present(defaultVal)) then
     nodeAs1dString = defaultVal
   else
@@ -1015,11 +1038,11 @@ function output_as1dString(self)  result(output)                   !ToDo: SR: Re
   class(tNode), pointer :: output_list
   integer :: o
 
-  output_list   => self%get('output',defaultVal=emptyList)
+  output_list => self%get('output',defaultVal=emptyList)
   allocate(output(output_list%length))
   do o = 1, output_list%length
     output(o) = output_list%get_asString(o)
-  enddo
+  end do
 
 end function output_as1dString
 
@@ -1042,7 +1065,7 @@ function tNode_get_byKey_asIndex(self,key)  result(keyIndex)
   do while (associated(item%next) .and. item%key /= key)
     item => item%next
     keyIndex = keyIndex+1
-  enddo
+  end do
 
   if (item%key /= key) call IO_error(140,ext_msg=key)
 
@@ -1087,7 +1110,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
     if (i /= 1) str = str//repeat(' ',indent_)
     str = str//'- '//item%node%asFormattedString(indent_+2)
     item => item%next
-  enddo
+  end do
 
 end function tList_asFormattedString
 
@@ -1119,9 +1142,9 @@ recursive function tDict_asFormattedString(self,indent) result(str)
         str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2)
       class default
         str = str//trim(item%key)//':'//IO_EOL//repeat(' ',indent_+2)//item%node%asFormattedString(indent_+2)
-    endselect
+    end select
     item => item%next
-  enddo
+  end do
 
 end function tDict_asFormattedString
 
@@ -1190,13 +1213,14 @@ function tList_as1dFloat(self)
   type(tItem),   pointer :: item
   type(tScalar), pointer :: scalar
 
+
   allocate(tList_as1dFloat(self%length))
   item => self%first
   do i = 1, self%length
     scalar => item%node%asScalar()
     tList_as1dFloat(i) = scalar%asFloat()
     item => item%next
-  enddo
+  end do
 
 end function tList_as1dFloat
 
@@ -1213,6 +1237,7 @@ function tList_as2dFloat(self)
   class(tNode), pointer :: row
   type(tList),  pointer :: row_data
 
+
   row => self%get(1)
   row_data => row%asList()
   allocate(tList_as2dFloat(self%length,row_data%length))
@@ -1220,9 +1245,9 @@ function tList_as2dFloat(self)
   do i=1,self%length
     row => self%get(i)
     row_data => row%asList()
-    if(row_data%length /= size(tList_as2dFloat,2)) call IO_error(709,ext_msg='Varying number of columns')
+    if (row_data%length /= size(tList_as2dFloat,2)) call IO_error(709,ext_msg='Varying number of columns')
     tList_as2dFloat(i,:) = self%get_as1dFloat(i)
-  enddo
+  end do
 
 end function tList_as2dFloat
 
@@ -1239,13 +1264,14 @@ function tList_as1dInt(self)
   type(tItem),   pointer :: item
   type(tScalar), pointer :: scalar
 
+
   allocate(tList_as1dInt(self%length))
   item => self%first
   do i = 1, self%length
     scalar => item%node%asScalar()
     tList_as1dInt(i) = scalar%asInt()
     item => item%next
-  enddo
+  end do
 
 end function tList_as1dInt
 
@@ -1262,13 +1288,14 @@ function tList_as1dBool(self)
   type(tItem),   pointer :: item
   type(tScalar), pointer :: scalar
 
+
   allocate(tList_as1dBool(self%length))
   item => self%first
   do i = 1, self%length
     scalar => item%node%asScalar()
     tList_as1dBool(i) = scalar%asBool()
     item => item%next
-  enddo
+  end do
 
 end function tList_as1dBool
 
@@ -1285,13 +1312,14 @@ function tList_as1dString(self)
   type(tItem),   pointer :: item
   type(tScalar), pointer :: scalar
 
+
   len_max = 0
   item => self%first
   do i = 1, self%length
     scalar => item%node%asScalar()
     len_max = max(len_max, len_trim(scalar%asString()))
     item => item%next
-  enddo
+  end do
 
   allocate(character(len=len_max) :: tList_as1dString(self%length))
   item => self%first
diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90
index 80cde388f..ee5fd4c82 100644
--- a/src/grid/spectral_utilities.f90
+++ b/src/grid/spectral_utilities.f90
@@ -156,7 +156,7 @@ subroutine spectral_utilities_init
   integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset
   integer(C_INTPTR_T), parameter :: &
     scalarSize = 1_C_INTPTR_T, &
-    vecSize    = 3_C_INTPTR_T, &
+    vectorSize = 3_C_INTPTR_T, &
     tensorSize = 9_C_INTPTR_T
   character(len=*), parameter :: &
     PETSCDEBUG = ' -snes_view -snes_monitor '
@@ -274,7 +274,7 @@ subroutine spectral_utilities_init
   call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, &
                    gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T ,              gridFFTW(2),local_K])     ! place a pointer for a fourier tensor representation
 
-  vectorField = fftw_alloc_complex(vecSize*alloc_local)
+  vectorField = fftw_alloc_complex(vectorSize*alloc_local)
   call c_f_pointer(vectorField, vectorField_real,   [3_C_INTPTR_T,&
                    2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K])     ! place a pointer for a real vector representation
   call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,&
@@ -288,42 +288,42 @@ subroutine spectral_utilities_init
 
 !--------------------------------------------------------------------------------------------------
 ! tensor MPI fftw plans
-  planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
-                                        tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
-                                                       tensorField_real, tensorField_fourier, &     ! input data, output data
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
+  planTensorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),tensorSize, &
+                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
+                                               tensorField_real,tensorField_fourier, &
+                                               PETSC_COMM_WORLD,FFTW_planner_flag)
   if (.not. c_associated(planTensorForth)) error stop 'FFTW error'
-  planTensorBack  = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
-                                       tensorSize,  FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
-                                                        tensorField_fourier,tensorField_real, &     ! input data, output data
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
+  planTensorBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),tensorSize, &
+                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
+                                               tensorField_fourier,tensorField_real, &
+                                               PETSC_COMM_WORLD, FFTW_planner_flag)
   if (.not. c_associated(planTensorBack)) error stop 'FFTW error'
 
 !--------------------------------------------------------------------------------------------------
 ! vector MPI fftw plans
-  planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
-                                           vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,&! no. of transforms, default iblock and oblock
-                                                       vectorField_real, vectorField_fourier, &     ! input data, output data
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
-  if (.not. C_ASSOCIATED(planVectorForth)) error stop 'FFTW error'
-  planVectorBack  = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)],  &         ! dimension, logical length in each dimension in reversed order
-                                         vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock
-                                                      vectorField_fourier,vectorField_real, &       ! input data, output data
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
-  if (.not. C_ASSOCIATED(planVectorBack)) error stop 'FFTW error'
+  planVectorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),vectorSize, &
+                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
+                                               vectorField_real,vectorField_fourier, &
+                                               PETSC_COMM_WORLD,FFTW_planner_flag)
+  if (.not. c_associated(planVectorForth)) error stop 'FFTW error'
+  planVectorBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),vectorSize, &
+                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
+                                               vectorField_fourier,vectorField_real, &
+                                               PETSC_COMM_WORLD, FFTW_planner_flag)
+  if (.not. c_associated(planVectorBack)) error stop 'FFTW error'
 
 !--------------------------------------------------------------------------------------------------
 ! scalar MPI fftw plans
-  planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
-                                       scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
-                                                       scalarField_real, scalarField_fourier, &     ! input data, output data
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
-  if (.not. C_ASSOCIATED(planScalarForth)) error stop 'FFTW error'
-  planScalarBack  = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order, no. of transforms
-                                       scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
-                                                       scalarField_fourier,scalarField_real, &      ! input data, output data
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
-  if (.not. C_ASSOCIATED(planScalarBack)) error stop 'FFTW error'
+  planScalarForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),scalarSize, &
+                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
+                                               scalarField_real,scalarField_fourier, &
+                                               PETSC_COMM_WORLD,FFTW_planner_flag)
+  if (.not. c_associated(planScalarForth)) error stop 'FFTW error'
+  planScalarBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),scalarSize, &
+                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
+                                               scalarField_fourier,scalarField_real, &
+                                               PETSC_COMM_WORLD, FFTW_planner_flag)
+  if (.not. c_associated(planScalarBack)) error stop 'FFTW error'
 
 !--------------------------------------------------------------------------------------------------
 ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
diff --git a/src/lattice.f90 b/src/lattice.f90
index 3089aa30b..e18c71edf 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -587,8 +587,8 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Non-schmid projections for bcc with up to 6 coefficients
-! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17)
-! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1
+! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17)
+! https://doi.org/10.1016/j.actamat.2008.07.037, table 1
 !--------------------------------------------------------------------------------------------------
 function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
 
@@ -602,6 +602,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
   type(rotation)                                       :: R
   integer                                              :: i
 
+
   if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix'
 
   coordinateSystem  = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,'cI',0.0_pReal)
@@ -634,7 +635,9 @@ end function lattice_nonSchmidMatrix
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Slip-slip interaction matrix
-!> details only active slip systems are considered
+!> @details only active slip systems are considered
+!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (fcc: Tab S4-1, bcc: Tab S5-1)
+!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hex: Tab 3b)
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
 
@@ -646,6 +649,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
   integer, dimension(:),   allocatable :: NslipMax
   integer, dimension(:,:), allocatable :: interactionTypes
 
+
   integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: &
     FCC_INTERACTIONSLIPSLIP = reshape( [&
        1, 2, 2, 4, 7, 5, 3, 5, 5, 4, 6, 7,  10,11,10,11,12,13, & ! -----> acting (forest)
@@ -750,41 +754,113 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
   integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
     HEX_INTERACTIONSLIPSLIP = reshape( [&
     ! basal      prism      1. pyr<a>           1. pyr<c+a>                           2. pyr<c+a>
-       1, 2, 2,   3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13,13,13,13,13,13,13,  21,21,21,21,21,21, & ! -----> acting (forest)
-       2, 1, 2,   3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13,13,13,13,13,13,13,  21,21,21,21,21,21, & ! | basal
-       2, 2, 1,   3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13,13,13,13,13,13,13,  21,21,21,21,21,21, & ! |
+       1, 2, 2,   3, 4, 4,   9,10, 9, 9,10, 9,  20,21,22,22,21,20,20,21,22,22,21,20,  47,47,48,47,47,48, & ! -----> acting (forest)
+       2, 1, 2,   4, 3, 4,  10, 9, 9,10, 9, 9,  22,22,21,20,20,21,22,22,21,20,20,21,  47,48,47,47,48,47, & ! | basal
+       2, 2, 1,   4, 4, 3,   9, 9,10, 9, 9,10,  21,20,20,21,22,22,21,20,20,21,22,22,  48,47,47,48,47,47, & ! |
                                                                                                            ! v
-       6, 6, 6,   4, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14,14,14,14,14,14,14,  22,22,22,22,22,22, & ! reacting (primary)
-       6, 6, 6,   5, 4, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14,14,14,14,14,14,14,  22,22,22,22,22,22, & ! prism
-       6, 6, 6,   5, 5, 4,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14,14,14,14,14,14,14,  22,22,22,22,22,22, &
+       7, 8, 8,   5, 6, 6,  11,12,11,11,12,11,  23,24,25,25,24,23,23,24,25,25,24,23,  49,49,50,49,49,50, & ! reacting (primary)
+       8, 7, 8,   6, 5, 6,  12,11,11,12,11,11,  25,25,24,23,23,24,25,25,24,23,23,24,  49,50,49,49,50,49, & ! prism
+       8, 8, 7,   6, 6, 5,  11,11,12,11,11,12,  24,23,23,24,25,25,24,23,23,24,25,25,  50,49,49,50,49,49, &
 
-      12,12,12,  11,11,11,   9,10,10,10,10,10,  15,15,15,15,15,15,15,15,15,15,15,15,  23,23,23,23,23,23, &
-      12,12,12,  11,11,11,  10, 9,10,10,10,10,  15,15,15,15,15,15,15,15,15,15,15,15,  23,23,23,23,23,23, &
-      12,12,12,  11,11,11,  10,10, 9,10,10,10,  15,15,15,15,15,15,15,15,15,15,15,15,  23,23,23,23,23,23, &
-      12,12,12,  11,11,11,  10,10,10, 9,10,10,  15,15,15,15,15,15,15,15,15,15,15,15,  23,23,23,23,23,23, & ! 1. pyr<a>
-      12,12,12,  11,11,11,  10,10,10,10, 9,10,  15,15,15,15,15,15,15,15,15,15,15,15,  23,23,23,23,23,23, &
-      12,12,12,  11,11,11,  10,10,10,10,10, 9,  15,15,15,15,15,15,15,15,15,15,15,15,  23,23,23,23,23,23, &
+      18,19,18,  16,17,16,  13,14,14,15,14,14,  26,26,27,28,28,27,29,29,27,28,28,27,  51,52,51,51,52,51, &
+      19,18,18,  17,16,16,  14,13,14,14,15,14,  28,27,26,26,27,28,28,27,29,29,27,28,  51,51,52,51,51,52, &
+      18,18,19,  16,16,17,  14,14,13,14,14,15,  27,28,28,27,26,26,27,28,28,27,29,29,  52,51,51,52,51,51, &
+      18,19,18,  16,17,16,  15,14,14,13,14,14,  29,29,27,28,28,27,26,26,27,28,28,27,  51,52,51,51,52,51, & ! 1. pyr<a>
+      19,18,18,  17,16,16,  14,15,14,14,13,14,  28,27,29,29,27,28,28,27,26,26,27,28,  51,51,52,51,51,52, &
+      18,18,19,  16,16,17,  14,14,15,14,14,13,  27,28,28,27,29,29,27,28,28,27,26,26,  52,51,51,52,51,51, &
 
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  16,17,17,17,17,17,17,17,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,16,17,17,17,17,17,17,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,16,17,17,17,17,17,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,16,17,17,17,17,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,16,17,17,17,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,16,17,17,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,17,16,17,17,17,17,17,  24,24,24,24,24,24, & ! 1. pyr<c+a>
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,17,17,16,17,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,17,17,17,16,17,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,17,17,17,17,16,17,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,17,17,17,17,17,16,17,  24,24,24,24,24,24, &
-      20,20,20,  19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,17,17,17,17,17,17,16,  24,24,24,24,24,24, &
+      44,45,46,  41,42,43,  37,38,39,40,38,39,  30,31,32,32,32,33,34,35,32,32,32,36,  53,54,55,53,54,56, &
+      46,45,44,  43,42,41,  37,39,38,40,39,38,  31,30,36,32,32,32,35,34,33,32,32,32,  56,54,53,55,54,53, &
+      45,46,44,  42,43,41,  39,37,38,39,40,38,  32,36,30,31,32,32,32,33,34,35,32,32,  56,53,54,55,53,54, &
+      45,44,46,  42,41,43,  38,37,39,38,40,39,  32,32,31,30,36,32,32,32,35,34,33,32,  53,56,54,53,55,54, &
+      46,44,45,  43,41,42,  38,39,37,38,39,40,  32,32,32,36,30,31,32,32,32,33,34,35,  54,56,53,54,55,53, &
+      44,46,45,  41,43,42,  39,38,37,39,38,40,  33,32,32,32,31,30,36,32,32,32,35,34,  54,53,56,54,53,55, &
+      44,45,46,  41,42,43,  40,38,39,37,38,39,  34,35,32,32,32,36,30,31,32,32,32,33,  53,54,56,53,54,55, & ! 1. pyr<c+a>
+      46,45,44,  43,42,41,  40,39,38,37,39,38,  35,34,33,32,32,32,31,30,36,32,32,32,  55,54,53,56,54,53, &
+      45,46,44,  42,43,41,  39,40,38,39,37,38,  32,33,34,35,32,32,32,36,30,31,32,32,  55,53,54,56,53,54, &
+      45,44,46,  42,41,43,  38,40,39,38,37,39,  32,32,35,34,33,32,32,32,31,30,36,32,  53,55,54,53,56,54, &
+      46,44,45,  43,41,42,  38,39,40,38,39,37,  32,32,32,33,34,35,32,32,32,36,30,31,  54,55,53,54,56,53, &
+      44,46,45,  41,43,42,  39,38,40,39,38,37,  36,32,32,32,35,34,33,32,32,32,31,30,  54,53,55,54,53,56, &
 
-      30,30,30,  29,29,29,  28,28,28,28,28,28,  27,27,27,27,27,27,27,27,27,27,27,27,  25,26,26,26,26,26, &
-      30,30,30,  29,29,29,  28,28,28,28,28,28,  27,27,27,27,27,27,27,27,27,27,27,27,  26,25,26,26,26,26, &
-      30,30,30,  29,29,29,  28,28,28,28,28,28,  27,27,27,27,27,27,27,27,27,27,27,27,  26,26,25,26,26,26, &
-      30,30,30,  29,29,29,  28,28,28,28,28,28,  27,27,27,27,27,27,27,27,27,27,27,27,  26,26,26,25,26,26, & ! 2. pyr<c+a>
-      30,30,30,  29,29,29,  28,28,28,28,28,28,  27,27,27,27,27,27,27,27,27,27,27,27,  26,26,26,26,25,26, &
-      30,30,30,  29,29,29,  28,28,28,28,28,28,  27,27,27,27,27,27,27,27,27,27,27,27,  26,26,26,26,26,25  &
+      68,68,69,  66,66,67,  64,64,65,64,65,65,  60,61,61,60,62,62,60,63,63,60,62,62,  57,58,58,59,58,58, &
+      68,69,68,  66,67,66,  65,64,64,65,64,64,  62,62,60,61,61,60,62,62,60,63,63,60,  58,57,58,58,59,58, &
+      69,68,68,  67,66,66,  64,65,64,64,65,64,  63,60,62,62,60,61,61,60,62,62,60,63,  58,58,57,58,58,59, &
+      68,68,69,  66,66,67,  64,64,65,64,64,65,  60,63,63,60,62,62,60,61,61,60,62,62,  59,58,58,57,58,58, & ! 2. pyr<c+a>
+      68,69,68,  66,67,66,  65,64,64,65,64,64,  62,62,60,63,63,60,62,62,60,61,61,60,  58,59,58,58,57,58, &
+      69,68,68,  67,66,66,  64,65,64,64,65,64,  61,60,62,62,60,63,63,60,62,62,60,61,  58,58,59,58,58,57  &
       ],shape(HEX_INTERACTIONSLIPSLIP))                                                             !< Slip-slip interaction types for hex (onion peel naming scheme)
+                                                                                                    !< 10.1016/j.ijplas.2014.06.010 table 3
+                                                                                                    !< 10.1080/14786435.2012.699689 table 2 and 3
+                                                                                                    !< index & label & description
+                                                                                                    !<  1 & S1 & basal self-interaction
+                                                                                                    !<  2 &  1 & basal/basal coplanar
+                                                                                                    !<  3 &  3 & basal/prismatic collinear
+                                                                                                    !<  4 &  4 & basal/prismatic non-collinear
+                                                                                                    !<  5 & S2 & prismatic self-interaction
+                                                                                                    !<  6 &  2 & prismatic/prismatic
+                                                                                                    !<  7 &  5 & prismatic/basal collinear
+                                                                                                    !<  8 &  6 & prismatic/basal non-collinear
+                                                                                                    !<  9 &  - & basal/pyramidal <a> non-collinear
+                                                                                                    !< 10 &  - & basal/pyramidal <a> collinear
+                                                                                                    !< 11 &  - & prismatic/pyramidal <a> non-collinear
+                                                                                                    !< 12 &  - & prismatic/pyramidal <a> collinear
+                                                                                                    !< 13 &  - & pyramidal <a> self-interaction
+                                                                                                    !< 14 &  - & pyramidal <a> non-collinear
+                                                                                                    !< 15 &  - & pyramidal <a> collinear
+                                                                                                    !< 16 &  - & pyramidal <a>/prismatic non-collinear
+                                                                                                    !< 17 &  - & pyramidal <a>/prismatic collinear
+                                                                                                    !< 18 &  - & pyramidal <a>/basal non-collinear
+                                                                                                    !< 19 &  - & pyramidal <a>/basal collinear
+                                                                                                    !< 20 &  - & basal/1. order pyramidal <c+a> semi-collinear
+                                                                                                    !< 21 &  - & basal/1. order pyramidal <c+a>
+                                                                                                    !< 22 &  - & basal/1. order pyramidal <c+a>
+                                                                                                    !< 23 &  - & prismatic/1. order pyramidal <c+a> semi-collinear
+                                                                                                    !< 24 &  - & prismatic/1. order pyramidal <c+a>
+                                                                                                    !< 25 &  - & prismatic/1. order pyramidal <c+a> semi-coplanar?
+                                                                                                    !< 26 &  - & pyramidal <a>/1. order pyramidal <c+a> coplanar
+                                                                                                    !< 27 &  - & pyramidal <a>/1. order pyramidal <c+a>
+                                                                                                    !< 28 &  - & pyramidal <a>/1. order pyramidal <c+a> semi-collinear
+                                                                                                    !< 29 &  - & pyramidal <a>/1. order pyramidal <c+a> semi-coplanar
+                                                                                                    !< 30 &  - & 1. order pyramidal <c+a> self-interaction
+                                                                                                    !< 31 &  - & 1. order pyramidal <c+a> coplanar
+                                                                                                    !< 32 &  - & 1. order pyramidal <c+a>
+                                                                                                    !< 33 &  - & 1. order pyramidal <c+a>
+                                                                                                    !< 34 &  - & 1. order pyramidal <c+a> semi-coplanar
+                                                                                                    !< 35 &  - & 1. order pyramidal <c+a> semi-coplanar
+                                                                                                    !< 36 &  - & 1. order pyramidal <c+a> collinear
+                                                                                                    !< 37 &  - & 1. order pyramidal <c+a>/pyramidal <a> coplanar
+                                                                                                    !< 38 &  - & 1. order pyramidal <c+a>/pyramidal <a> semi-collinear
+                                                                                                    !< 39 &  - & 1. order pyramidal <c+a>/pyramidal <a>
+                                                                                                    !< 40 &  - & 1. order pyramidal <c+a>/pyramidal <a> semi-coplanar
+                                                                                                    !< 41 &  - & 1. order pyramidal <c+a>/prismatic semi-collinear
+                                                                                                    !< 42 &  - & 1. order pyramidal <c+a>/prismatic semi-coplanar
+                                                                                                    !< 43 &  - & 1. order pyramidal <c+a>/prismatic
+                                                                                                    !< 44 &  - & 1. order pyramidal <c+a>/basal semi-collinear
+                                                                                                    !< 45 &  - & 1. order pyramidal <c+a>/basal
+                                                                                                    !< 46 &  - & 1. order pyramidal <c+a>/basal
+                                                                                                    !< 47 &  8 & basal/2. order pyramidal <c+a> non-collinear
+                                                                                                    !< 48 &  7 & basal/2. order pyramidal <c+a> semi-collinear
+                                                                                                    !< 49 & 10 & prismatic/2. order pyramidal <c+a>
+                                                                                                    !< 50 &  9 & prismatic/2. order pyramidal <c+a> semi-collinear
+                                                                                                    !< 51 &  - & pyramidal <a>/2. order pyramidal <c+a>
+                                                                                                    !< 52 &  - & pyramidal <a>/2. order pyramidal <c+a> semi collinear
+                                                                                                    !< 53 &  - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
+                                                                                                    !< 54 &  - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
+                                                                                                    !< 55 &  - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
+                                                                                                    !< 56 &  - & 1. order pyramidal <c+a>/2. order pyramidal <c+a> collinear
+                                                                                                    !< 57 & S3 & 2. order pyramidal <c+a> self-interaction
+                                                                                                    !< 58 & 16 & 2. order pyramidal <c+a> non-collinear
+                                                                                                    !< 59 & 15 & 2. order pyramidal <c+a> semi-collinear
+                                                                                                    !< 60 &  - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
+                                                                                                    !< 61 &  - & 2. order pyramidal <c+a>/1. order pyramidal <c+a> collinear
+                                                                                                    !< 62 &  - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
+                                                                                                    !< 63 &  - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
+                                                                                                    !< 64 &  - & 2. order pyramidal <c+a>/pyramidal <a> non-collinear
+                                                                                                    !< 65 &  - & 2. order pyramidal <c+a>/pyramidal <a> semi-collinear
+                                                                                                    !< 66 & 14 & 2. order pyramidal <c+a>/prismatic non-collinear
+                                                                                                    !< 67 & 13 & 2. order pyramidal <c+a>/prismatic semi-collinear
+                                                                                                    !< 68 & 12 & 2. order pyramidal <c+a>/basal non-collinear
+                                                                                                    !< 69 & 11 & 2. order pyramidal <c+a>/basal semi-collinear
 
   integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: &
     BCT_INTERACTIONSLIPSLIP = reshape( [&
diff --git a/src/material.f90 b/src/material.f90
index 52940ee4a..1e3a4b4ec 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -108,8 +108,14 @@ subroutine parse()
   homogenizations => config_material%get('homogenization')
 
   call sanityCheck(materials, homogenizations)
+
+#if defined (__GFORTRAN__)
   material_name_phase          = getKeys(phases)
   material_name_homogenization = getKeys(homogenizations)
+#else
+  material_name_phase          = phases%Keys()
+  material_name_homogenization = homogenizations%Keys()
+#endif
 
   allocate(homogenization_Nconstituents(homogenizations%length))
   do h=1, homogenizations%length
@@ -203,9 +209,9 @@ subroutine sanityCheck(materials,homogenizations)
 
 end subroutine sanityCheck
 
-
+#if defined (__GFORTRAN__)
 !--------------------------------------------------------------------------------------------------
-!> @brief Get all keys from a dictionary
+!> @brief %keys() is broken on gfortran
 !--------------------------------------------------------------------------------------------------
 function getKeys(dict)
 
@@ -228,5 +234,6 @@ function getKeys(dict)
   end do
 
 end function getKeys
+#endif
 
 end module material
diff --git a/src/math.f90 b/src/math.f90
index 3ea4d6a08..db0666ab0 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -84,38 +84,34 @@ contains
 subroutine math_init
 
   real(pReal), dimension(4) :: randTest
-  integer :: &
-    randSize, &
-    randomSeed                                                                                      !< fixed seeding for pseudo-random number generator, Default 0: use random seed
-  integer, dimension(:), allocatable :: randInit
+  integer :: randSize
+  integer, dimension(:), allocatable :: seed
   class(tNode), pointer :: &
     num_generic
 
+
   print'(/,1x,a)', '<<<+-  math init  -+>>>'; flush(IO_STDOUT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
-  randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)
 
   call random_seed(size=randSize)
-  allocate(randInit(randSize))
-  if (randomSeed > 0) then
-    randInit = randomSeed
+  allocate(seed(randSize))
+
+  if (num_generic%contains('random_seed')) then
+    seed = num_generic%get_as1dInt('random_seed',requiredSize=randSize)
   else
     call random_seed()
-    call random_seed(get = randInit)
-    randInit(2:randSize) = randInit(1)
-  endif
+    call random_seed(get = seed)
+  end if
 
-  call random_seed(put = randInit)
+  call random_seed(put = seed)
   call random_number(randTest)
 
-  print'(/,a,i2)',                ' size  of random seed:     ', randSize
-  print'(  a,i0)',                ' value of random seed:     ', randInit(1)
-  print'(  a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
+  print'(/,a,i2)',              ' size  of random seed:     ', randSize
+  print*,                       'value of random seed:     ', seed
+  print'(  a,4(/,26x,f17.14))', ' start of random sequence: ', randTest
 
-  call random_seed(put = randInit)
-
-  call selfTest
+  call selfTest()
 
 end subroutine math_init
 
diff --git a/src/rotations.f90 b/src/rotations.f90
index b4728e38b..b73fbe8da 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -270,7 +270,7 @@ pure elemental subroutine standardize(self)
 
   class(rotation), intent(inout) :: self
 
-  if (self%q(1) < 0.0_pReal) self%q = - self%q
+  if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
 
 end subroutine standardize
 
@@ -450,7 +450,7 @@ pure function qu2om(qu) result(om)
   om(3,2) = 2.0_pReal*(qu(4)*qu(3)+qu(1)*qu(2))
   om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
 
-  if (P < 0.0_pReal) om = transpose(om)
+  if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
 
 end function qu2om
 
@@ -480,7 +480,7 @@ pure function qu2eu(qu) result(eu)
           atan2( 2.0_pReal*chi, q03-q12 ), &
           atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
   endif degenerated
-  where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
+  where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
 
 end function qu2eu
 
@@ -602,7 +602,7 @@ pure function om2qu(om) result(qu)
           qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
       endif
   endif
-  if(qu(1)<0._pReal) qu =-1.0_pReal * qu
+  if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
   qu = qu*[1.0_pReal,P,P,P]
 
 end function om2qu
@@ -628,7 +628,7 @@ pure function om2eu(om) result(eu)
     eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
   end if
   where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
-  where(eu<0.0_pReal)          eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
+  where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
 
 end function om2eu
 
@@ -735,7 +735,7 @@ pure function eu2qu(eu) result(qu)
         -P*sPhi*cos(ee(1)-ee(3)), &
         -P*sPhi*sin(ee(1)-ee(3)), &
         -P*cPhi*sin(ee(1)+ee(3))]
-  if(qu(1) < 0.0_pReal) qu = qu * (-1.0_pReal)
+  if(sign(1.0_pReal,qu(1)) < 0.0_pReal) qu = qu * (-1.0_pReal)
 
 end function eu2qu
 
@@ -792,7 +792,7 @@ pure function eu2ax(eu) result(ax)
   else
     ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ]                                   ! passive axis-angle pair so a minus sign in front
     ax(4) = alpha
-    if (alpha < 0.0_pReal) ax = -ax                                                                 ! ensure alpha is positive
+    if (sign(1.0_pReal,alpha) < 0.0_pReal) ax = -ax                                                 ! ensure alpha is positive
   end if
 
 end function eu2ax