diff --git a/PRIVATE b/PRIVATE
index 1766c855e..0e82975b2 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 1766c855e54668d6225c9b22db536be7052605bc
+Subproject commit 0e82975b23a1bd4310c523388f1cadf1b8e03dd0
diff --git a/VERSION b/VERSION
index 435d3544f..00e433ed7 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-3.0.0-alpha6-371-g18d862cdb
+3.0.0-alpha6-395-g9696c1967
diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py
index 725567b7f..7a06a2c69 100644
--- a/python/damask/_configmaterial.py
+++ b/python/damask/_configmaterial.py
@@ -482,7 +482,7 @@ class ConfigMaterial(Config):
         """
         N,n,shaped = 1,1,{}
 
-        map_dim = {'O':-1,'F_i':-2}
+        map_dim = {'O':-1,'V_e':-2}
         for k,v in kwargs.items():
             shaped[k] = np.array(v)
             s = shaped[k].shape[:map_dim.get(k,None)]
@@ -494,12 +494,12 @@ class ConfigMaterial(Config):
         if 'v' not in kwargs:
             shaped['v'] = np.broadcast_to(1/n,(N,n))
 
-        map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)}
+        map_shape = {'O':(N,n,4),'V_e':(N,n,3,3)}
         for k,v in shaped.items():
             target = map_shape.get(k,(N,n))
             obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape, target, mode = 'right')), target)
             for i in range(N):
-                if k in ['phase','O','v','F_i']:
+                if k in ['phase','O','v','V_e']:
                     for j in range(n):
                         mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j]
                 else:
diff --git a/python/damask/_grid.py b/python/damask/_grid.py
index f7acfe5aa..4e82a4180 100644
--- a/python/damask/_grid.py
+++ b/python/damask/_grid.py
@@ -50,15 +50,16 @@ class Grid:
             Coordinates of grid origin in meter. Defaults to [0.0,0.0,0.0].
         initial_conditions : dictionary, optional
             Labels and values of the inital conditions at each material point.
-        comments : (list of) str, optional
-            Comments, e.g. history of operations.
+        comments : str or iterable of str, optional
+            Additional, human-readable information, e.g. history of operations.
 
         """
         self.material = material
         self.size = size                                                        # type: ignore
         self.origin = origin                                                    # type: ignore
         self.initial_conditions = {} if initial_conditions is None else initial_conditions
-        self.comments = [] if comments is None else comments                    # type: ignore
+        comments_ = [comments] if isinstance(comments,str) else comments
+        self.comments = [] if comments_ is None else [str(c) for c in comments_]
 
     def __repr__(self) -> str:
         """Give short human-readable summary."""
diff --git a/python/damask/_result.py b/python/damask/_result.py
index bf58ec160..1ad25bdbf 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -1,4 +1,5 @@
 import multiprocessing as mp
+from multiprocessing.synchronize import Lock
 import re
 import fnmatch
 import os
@@ -10,6 +11,7 @@ from pathlib import Path
 from functools import partial
 from collections import defaultdict
 from collections.abc import Iterable
+from typing import Union, Callable, Any, Sequence, Literal, Dict, List, Tuple
 
 import h5py
 import numpy as np
@@ -22,19 +24,21 @@ from . import grid_filters
 from . import mechanics
 from . import tensor
 from . import util
+from ._typehints import FloatSequence, IntSequence
 
 h5py3 = h5py.__version__[0] == '3'
 
 chunk_size = 1024**2//8                                                                             # for compression in HDF5
 prefix_inc = 'increment_'
 
-def _read(dataset):
+def _read(dataset: h5py._hl.dataset.Dataset) -> np.ndarray:
     """Read a dataset and its metadata into a numpy.ndarray."""
     metadata = {k:(v.decode() if not h5py3 and type(v) is bytes else v) for k,v in dataset.attrs.items()}
-    dtype = np.dtype(dataset.dtype,metadata=metadata)
+    dtype = np.dtype(dataset.dtype,metadata=metadata)                                               # type: ignore
     return np.array(dataset,dtype=dtype)
 
-def _match(requested,existing):
+def _match(requested,
+           existing: h5py._hl.base.KeysViewHDF5) -> List[Any]:
     """Find matches among two sets of labels."""
     def flatten_list(list_of_lists):
         return [e for e_ in list_of_lists for e in e_]
@@ -50,7 +54,10 @@ def _match(requested,existing):
     return sorted(set(flatten_list([fnmatch.filter(existing,r) for r in requested_])),
                   key=util.natural_sort)
 
-def _empty_like(dataset,N_materialpoints,fill_float,fill_int):
+def _empty_like(dataset: np.ma.core.MaskedArray,
+                N_materialpoints: int,
+                fill_float: float,
+                fill_int: int) -> np.ma.core.MaskedArray:
     """Create empty numpy.ma.MaskedArray."""
     return ma.array(np.empty((N_materialpoints,)+dataset.shape[1:],dataset.dtype),
                     fill_value = fill_float if dataset.dtype in np.sctypes['float'] else fill_int,
@@ -84,7 +91,7 @@ class Result:
 
     """
 
-    def __init__(self,fname):
+    def __init__(self, fname: Union[str, Path]):
         """
         New result view bound to a HDF5 file.
 
@@ -102,7 +109,7 @@ class Result:
             if self.version_major != 0 or not 12 <= self.version_minor <= 14:
                 raise TypeError(f'unsupported DADF5 version "{self.version_major}.{self.version_minor}"')
             if self.version_major == 0 and self.version_minor < 14:
-                self.export_setup = None
+                self.export_setup = None                                                            # type: ignore
 
             self.structured = 'cells' in f['geometry'].attrs.keys()
 
@@ -111,7 +118,7 @@ class Result:
                 self.size   = f['geometry'].attrs['size']
                 self.origin = f['geometry'].attrs['origin']
             else:
-                self.add_curl = self.add_divergence = self.add_gradient = None
+                self.add_curl = self.add_divergence = self.add_gradient = None                      # type: ignore
 
             r = re.compile(rf'{prefix_inc}([0-9]+)')
             self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort)
@@ -126,7 +133,7 @@ class Result:
             self.phase           = f['cell_to/phase']['label'].astype('str')
             self.phases          = sorted(np.unique(self.phase),key=util.natural_sort)
 
-            self.fields = []
+            self.fields: List[str] = []
             for c in self.phases:
                 self.fields += f['/'.join([self.increments[0],'phase',c])].keys()
             for m in self.homogenizations:
@@ -144,14 +151,14 @@ class Result:
         self._protected = True
 
 
-    def __copy__(self):
+    def __copy__(self) -> "Result":
         """Create deep copy."""
         return copy.deepcopy(self)
 
     copy = __copy__
 
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """Give short human-readable summary."""
         with h5py.File(self.fname,'r') as f:
             header = [f'Created by {f.attrs["creator"]}',
@@ -171,12 +178,12 @@ class Result:
 
 
     def _manage_view(self,
-                     action,
-                     increments=None,
-                     times=None,
-                     phases=None,
-                     homogenizations=None,
-                     fields=None):
+                     action: Literal['set', 'add', 'del'],
+                     increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
+                     times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
+                     phases: Union[str, Sequence[str], bool] = None,
+                     homogenizations: Union[str, Sequence[str], bool] = None,
+                     fields: Union[str, Sequence[str], bool] = None) -> "Result":
         """
         Manages the visibility of the groups.
 
@@ -204,8 +211,7 @@ class Result:
                 datasets = '*'
             elif datasets is False:
                 datasets = []
-            choice = list(datasets).copy() if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \
-                    [datasets]
+            choice = [datasets] if not hasattr(datasets,'__iter__') or isinstance(datasets,str) else list(datasets) # type: ignore
 
             if   what == 'increments':
                 choice = [c if isinstance(c,str) and c.startswith(prefix_inc) else
@@ -216,7 +222,7 @@ class Result:
                 if choice == ['*']:
                     choice = self.increments
                 else:
-                    iterator = map(float,choice)
+                    iterator = map(float,choice)                                                    # type: ignore
                     choice = []
                     for c in iterator:
                         idx = np.searchsorted(self.times,c)
@@ -224,7 +230,7 @@ class Result:
                         if   np.isclose(c,self.times[idx]):
                             choice.append(self.increments[idx])
                         elif np.isclose(c,self.times[idx+1]):
-                            choice.append(self.increments[idx+1])
+                            choice.append(self.increments[idx+1])                                   # type: ignore
 
             valid = _match(choice,getattr(self,what))
             existing = set(self.visible[what])
@@ -241,7 +247,9 @@ class Result:
         return dup
 
 
-    def increments_in_range(self,start=None,end=None):
+    def increments_in_range(self,
+                            start: Union[str, int] = None,
+                            end: Union[str, int] = None) -> Sequence[int]:
         """
         Get all increments within a given range.
 
@@ -263,7 +271,9 @@ class Result:
                    self.incs[-1] if  end  is None else  end))
         return [i for i in self.incs if s <= i <= e]
 
-    def times_in_range(self,start=None,end=None):
+    def times_in_range(self,
+                       start: float = None,
+                       end: float = None) -> Sequence[int]:
         """
         Get all increments within a given time range.
 
@@ -286,12 +296,12 @@ class Result:
 
 
     def view(self,*,
-             increments=None,
-             times=None,
-             phases=None,
-             homogenizations=None,
-             fields=None,
-             protected=None):
+             increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
+             times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
+             phases: Union[str, Sequence[str], bool] = None,
+             homogenizations: Union[str, Sequence[str], bool] = None,
+             fields: Union[str, Sequence[str], bool] = None,
+             protected: bool = None) -> "Result":
         """
         Set view.
 
@@ -343,11 +353,11 @@ class Result:
 
 
     def view_more(self,*,
-                  increments=None,
-                  times=None,
-                  phases=None,
-                  homogenizations=None,
-                  fields=None):
+                  increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
+                  times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
+                  phases: Union[str, Sequence[str], bool] = None,
+                  homogenizations: Union[str, Sequence[str], bool] = None,
+                  fields: Union[str, Sequence[str], bool] = None) -> "Result":
         """
         Add to view.
 
@@ -386,11 +396,11 @@ class Result:
 
 
     def view_less(self,*,
-                  increments=None,
-                  times=None,
-                  phases=None,
-                  homogenizations=None,
-                  fields=None):
+                  increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
+                  times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
+                  phases: Union[str, Sequence[str], bool] = None,
+                  homogenizations: Union[str, Sequence[str], bool] = None,
+                  fields: Union[str, Sequence[str], bool] = None) -> "Result":
         """
         Remove from view.
 
@@ -427,7 +437,9 @@ class Result:
         return self._manage_view('del',increments,times,phases,homogenizations,fields)
 
 
-    def rename(self,name_src,name_dst):
+    def rename(self,
+               name_src: str,
+               name_dst: str):
         """
         Rename/move datasets (within the same group/folder).
 
@@ -468,7 +480,7 @@ class Result:
                                 del f[path_src]
 
 
-    def remove(self,name):
+    def remove(self, name: str):
         """
         Remove/delete datasets.
 
@@ -502,7 +514,7 @@ class Result:
                             if path in f.keys(): del f[path]
 
 
-    def list_data(self):
+    def list_data(self) -> List[str]:
         """
         Collect information on all active datasets in the file.
 
@@ -533,7 +545,8 @@ class Result:
         return msg
 
 
-    def enable_user_function(self,func):
+    def enable_user_function(self,
+                             func: Callable):
         globals()[func.__name__]=func
         print(f'Function {func.__name__} enabled in add_calculation.')
 
@@ -544,7 +557,7 @@ class Result:
 
 
     @property
-    def coordinates0_point(self):
+    def coordinates0_point(self) -> np.ndarray:
         """Initial/undeformed cell center coordinates."""
         if self.structured:
             return grid_filters.coordinates0_point(self.cells,self.size,self.origin).reshape(-1,3,order='F')
@@ -553,7 +566,7 @@ class Result:
                 return f['geometry/x_p'][()]
 
     @property
-    def coordinates0_node(self):
+    def coordinates0_node(self) -> np.ndarray:
         """Initial/undeformed nodal coordinates."""
         if self.structured:
             return grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
@@ -562,7 +575,7 @@ class Result:
                 return f['geometry/x_n'][()]
 
     @property
-    def geometry0(self):
+    def geometry0(self) -> VTK:
         """Initial/undeformed geometry."""
         if self.structured:
             return VTK.from_image_data(self.cells,self.size,self.origin)
@@ -575,7 +588,7 @@ class Result:
 
 
     @staticmethod
-    def _add_absolute(x):
+    def _add_absolute(x: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  np.abs(x['data']),
                 'label': f'|{x["label"]}|',
@@ -585,7 +598,7 @@ class Result:
                           'creator':     'add_absolute'
                           }
                  }
-    def add_absolute(self,x):
+    def add_absolute(self, x: str):
         """
         Add absolute value.
 
@@ -599,7 +612,7 @@ class Result:
 
 
     @staticmethod
-    def _add_calculation(**kwargs):
+    def _add_calculation(**kwargs) -> Dict[str, Any]:
         formula = kwargs['formula']
         for d in re.findall(r'#(.*?)#',formula):
             formula = formula.replace(f'#{d}#',f"kwargs['{d}']['data']")
@@ -617,7 +630,11 @@ class Result:
                           'creator':     'add_calculation'
                           }
                  }
-    def add_calculation(self,formula,name,unit='n/a',description=None):
+    def add_calculation(self,
+                        formula: str,
+                        name: str,
+                        unit: str = 'n/a',
+                        description: str = None):
         """
         Add result of a general formula.
 
@@ -661,13 +678,13 @@ class Result:
         ...                   'Mises equivalent of the Cauchy stress')
 
         """
-        dataset_mapping  = {d:d for d in set(re.findall(r'#(.*?)#',formula))}                       # datasets used in the formula
+        dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))}                        # datasets used in the formula
         args             = {'formula':formula,'label':name,'unit':unit,'description':description}
         self._add_generic_pointwise(self._add_calculation,dataset_mapping,args)
 
 
     @staticmethod
-    def _add_stress_Cauchy(P,F):
+    def _add_stress_Cauchy(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  mechanics.stress_Cauchy(P['data'],F['data']),
                 'label': 'sigma',
@@ -679,7 +696,9 @@ class Result:
                           'creator':     'add_stress_Cauchy'
                           }
                 }
-    def add_stress_Cauchy(self,P='P',F='F'):
+    def add_stress_Cauchy(self,
+                          P: str = 'P',
+                          F: str = 'F'):
         """
         Add Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
 
@@ -695,7 +714,7 @@ class Result:
 
 
     @staticmethod
-    def _add_determinant(T):
+    def _add_determinant(T: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  np.linalg.det(T['data']),
                 'label': f"det({T['label']})",
@@ -705,7 +724,7 @@ class Result:
                           'creator':     'add_determinant'
                           }
                 }
-    def add_determinant(self,T):
+    def add_determinant(self, T: str):
         """
         Add the determinant of a tensor.
 
@@ -727,7 +746,7 @@ class Result:
 
 
     @staticmethod
-    def _add_deviator(T):
+    def _add_deviator(T: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  tensor.deviatoric(T['data']),
                 'label': f"s_{T['label']}",
@@ -737,7 +756,7 @@ class Result:
                           'creator':     'add_deviator'
                           }
                  }
-    def add_deviator(self,T):
+    def add_deviator(self, T: str):
         """
         Add the deviatoric part of a tensor.
 
@@ -759,13 +778,15 @@ class Result:
 
 
     @staticmethod
-    def _add_eigenvalue(T_sym,eigenvalue):
+    def _add_eigenvalue(T_sym: Dict[str, Any], eigenvalue: Literal['max, mid, min']) -> Dict[str, Any]:
         if   eigenvalue == 'max':
             label,p = 'maximum',2
         elif eigenvalue == 'mid':
             label,p = 'intermediate',1
         elif eigenvalue == 'min':
             label,p = 'minimum',0
+        else:
+            raise ValueError(f'invalid eigenvalue: {eigenvalue}')
 
         return {
                 'data': tensor.eigenvalues(T_sym['data'])[:,p],
@@ -776,7 +797,9 @@ class Result:
                           'creator':     'add_eigenvalue'
                          }
                 }
-    def add_eigenvalue(self,T_sym,eigenvalue='max'):
+    def add_eigenvalue(self,
+                       T_sym: str,
+                       eigenvalue: Literal['max', 'mid', 'min'] = 'max'):
         """
         Add eigenvalues of symmetric tensor.
 
@@ -800,13 +823,16 @@ class Result:
 
 
     @staticmethod
-    def _add_eigenvector(T_sym,eigenvalue):
+    def _add_eigenvector(T_sym: Dict[str, Any], eigenvalue: Literal['max', 'mid', 'min']) -> Dict[str, Any]:
         if   eigenvalue == 'max':
             label,p = 'maximum',2
         elif eigenvalue == 'mid':
             label,p = 'intermediate',1
         elif eigenvalue == 'min':
             label,p = 'minimum',0
+        else:
+            raise ValueError(f'invalid eigenvalue: {eigenvalue}')
+
         return {
                 'data': tensor.eigenvectors(T_sym['data'])[:,p],
                 'label': f"v_{eigenvalue}({T_sym['label']})",
@@ -817,7 +843,9 @@ class Result:
                           'creator':     'add_eigenvector'
                          }
                }
-    def add_eigenvector(self,T_sym,eigenvalue='max'):
+    def add_eigenvector(self,
+                        T_sym: str,
+                        eigenvalue: Literal['max', 'mid', 'min'] = 'max'):
         """
         Add eigenvector of symmetric tensor.
 
@@ -834,7 +862,7 @@ class Result:
 
 
     @staticmethod
-    def _add_IPF_color(l,q):
+    def _add_IPF_color(l: FloatSequence, q: Dict[str, Any]) -> Dict[str, Any]:
         m = util.scale_to_coprime(np.array(l))
         lattice =  q['meta']['lattice']
         o = Orientation(rotation = q['data'],lattice=lattice)
@@ -849,7 +877,9 @@ class Result:
                           'creator':     'add_IPF_color'
                          }
                }
-    def add_IPF_color(self,l,q='O'):
+    def add_IPF_color(self,
+                      l: FloatSequence,
+                      q: str = 'O'):
         """
         Add RGB color tuple of inverse pole figure (IPF) color.
 
@@ -874,7 +904,7 @@ class Result:
 
 
     @staticmethod
-    def _add_maximum_shear(T_sym):
+    def _add_maximum_shear(T_sym: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  mechanics.maximum_shear(T_sym['data']),
                 'label': f"max_shear({T_sym['label']})",
@@ -884,7 +914,7 @@ class Result:
                           'creator':     'add_maximum_shear'
                           }
                  }
-    def add_maximum_shear(self,T_sym):
+    def add_maximum_shear(self, T_sym: str):
         """
         Add maximum shear components of symmetric tensor.
 
@@ -898,7 +928,7 @@ class Result:
 
 
     @staticmethod
-    def _add_equivalent_Mises(T_sym,kind):
+    def _add_equivalent_Mises(T_sym: Dict[str, Any], kind: str) -> Dict[str, Any]:
         k = kind
         if k is None:
             if T_sym['meta']['unit'] == '1':
@@ -918,7 +948,9 @@ class Result:
                           'creator':     'add_Mises'
                           }
                 }
-    def add_equivalent_Mises(self,T_sym,kind=None):
+    def add_equivalent_Mises(self,
+                             T_sym: str,
+                             kind: str = None):
         """
         Add the equivalent Mises stress or strain of a symmetric tensor.
 
@@ -949,10 +981,10 @@ class Result:
 
 
     @staticmethod
-    def _add_norm(x,ord):
+    def _add_norm(x: Dict[str, Any], ord: Union[int, float, Literal['fro', 'nuc']]) ->  Dict[str, Any]:
         o = ord
         if len(x['data'].shape) == 2:
-            axis = 1
+            axis: Union[int, Tuple[int, int]] = 1
             t = 'vector'
             if o is None: o = 2
         elif len(x['data'].shape) == 3:
@@ -971,7 +1003,9 @@ class Result:
                           'creator':     'add_norm'
                           }
                  }
-    def add_norm(self,x,ord=None):
+    def add_norm(self,
+                 x: str,
+                 ord: Union[int, float, Literal['fro', 'nuc']] = None):
         """
         Add the norm of vector or tensor.
 
@@ -987,7 +1021,7 @@ class Result:
 
 
     @staticmethod
-    def _add_stress_second_Piola_Kirchhoff(P,F):
+    def _add_stress_second_Piola_Kirchhoff(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']),
                 'label': 'S',
@@ -999,7 +1033,9 @@ class Result:
                           'creator':     'add_stress_second_Piola_Kirchhoff'
                           }
                 }
-    def add_stress_second_Piola_Kirchhoff(self,P='P',F='F'):
+    def add_stress_second_Piola_Kirchhoff(self,
+                                          P: str = 'P',
+                                          F: str = 'F'):
         """
         Add second Piola-Kirchhoff stress calculated from first Piola-Kirchhoff stress and deformation gradient.
 
@@ -1023,7 +1059,11 @@ class Result:
 
 
     @staticmethod
-    def _add_pole(q,uvw,hkl,with_symmetry,normalize):
+    def _add_pole(q: Dict[str, Any],
+                  uvw: FloatSequence,
+                  hkl: FloatSequence,
+                  with_symmetry: bool,
+                  normalize: bool) -> Dict[str, Any]:
         c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1
         brackets = ['[]','()','⟨⟩','{}'][(uvw is None)*1+with_symmetry*2]
         label = 'p^' + '{}{} {} {}{}'.format(brackets[0],
@@ -1042,7 +1082,13 @@ class Result:
                           'creator':     'add_pole'
                           }
                 }
-    def add_pole(self,q='O',*,uvw=None,hkl=None,with_symmetry=False,normalize=True):
+    def add_pole(self,
+                 q: str = 'O',
+                 *,
+                 uvw: FloatSequence = None,
+                 hkl: FloatSequence = None,
+                 with_symmetry: bool = False,
+                 normalize: bool = True):
         """
         Add lab frame vector along lattice direction [uvw] or plane normal (hkl).
 
@@ -1067,7 +1113,7 @@ class Result:
 
 
     @staticmethod
-    def _add_rotation(F):
+    def _add_rotation(F: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  mechanics.rotation(F['data']).as_matrix(),
                 'label': f"R({F['label']})",
@@ -1077,7 +1123,7 @@ class Result:
                           'creator':     'add_rotation'
                           }
                  }
-    def add_rotation(self,F):
+    def add_rotation(self, F: str):
         """
         Add rotational part of a deformation gradient.
 
@@ -1099,7 +1145,7 @@ class Result:
 
 
     @staticmethod
-    def _add_spherical(T):
+    def _add_spherical(T: Dict[str, Any]) -> Dict[str, Any]:
         return {
                 'data':  tensor.spherical(T['data'],False),
                 'label': f"p_{T['label']}",
@@ -1109,7 +1155,7 @@ class Result:
                           'creator':     'add_spherical'
                           }
                  }
-    def add_spherical(self,T):
+    def add_spherical(self, T: str):
         """
         Add the spherical (hydrostatic) part of a tensor.
 
@@ -1131,7 +1177,7 @@ class Result:
 
 
     @staticmethod
-    def _add_strain(F,t,m):
+    def _add_strain(F: Dict[str, Any], t: Literal['V', 'U'], m: float) -> Dict[str, Any]:
         return {
                 'data':  mechanics.strain(F['data'],t,m),
                 'label': f"epsilon_{t}^{m}({F['label']})",
@@ -1141,7 +1187,10 @@ class Result:
                           'creator':     'add_strain'
                           }
                  }
-    def add_strain(self,F='F',t='V',m=0.0):
+    def add_strain(self,
+                   F: str = 'F',
+                   t: Literal['V', 'U'] = 'V',
+                   m: float = 0.0):
         """
         Add strain tensor of a deformation gradient.
 
@@ -1177,7 +1226,7 @@ class Result:
 
 
     @staticmethod
-    def _add_stretch_tensor(F,t):
+    def _add_stretch_tensor(F: Dict[str, Any], t: str) -> Dict[str, Any]:
         return {
                 'data':  (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']),
                 'label': f"{t}({F['label']})",
@@ -1188,7 +1237,9 @@ class Result:
                           'creator':     'add_stretch_tensor'
                           }
                  }
-    def add_stretch_tensor(self,F='F',t='V'):
+    def add_stretch_tensor(self,
+                           F: str = 'F',
+                           t: Literal['V', 'U'] = 'V'):
         """
         Add stretch tensor of a deformation gradient.
 
@@ -1205,7 +1256,7 @@ class Result:
 
 
     @staticmethod
-    def _add_curl(f,size):
+    def _add_curl(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
         return {
                 'data':  grid_filters.curl(size,f['data']),
                 'label': f"curl({f['label']})",
@@ -1215,7 +1266,7 @@ class Result:
                           'creator':     'add_curl'
                           }
                  }
-    def add_curl(self,f):
+    def add_curl(self, f: str):
         """
         Add curl of a field.
 
@@ -1234,7 +1285,7 @@ class Result:
 
 
     @staticmethod
-    def _add_divergence(f,size):
+    def _add_divergence(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
         return {
                 'data':  grid_filters.divergence(size,f['data']),
                 'label': f"divergence({f['label']})",
@@ -1244,7 +1295,7 @@ class Result:
                           'creator':     'add_divergence'
                           }
                  }
-    def add_divergence(self,f):
+    def add_divergence(self, f: str):
         """
         Add divergence of a field.
 
@@ -1263,7 +1314,7 @@ class Result:
 
 
     @staticmethod
-    def _add_gradient(f,size):
+    def _add_gradient(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
         return {
                 'data':  grid_filters.gradient(size,f['data'] if len(f['data'].shape) == 4 else \
                                                     f['data'].reshape(f['data'].shape+(1,))),
@@ -1274,7 +1325,7 @@ class Result:
                           'creator':     'add_gradient'
                           }
                  }
-    def add_gradient(self,f):
+    def add_gradient(self, f: str):
         """
         Add gradient of a field.
 
@@ -1292,7 +1343,11 @@ class Result:
         self._add_generic_grid(self._add_gradient,{'f':f},{'size':self.size})
 
 
-    def _add_generic_grid(self,func,datasets,args={},constituents=None):
+    def _add_generic_grid(self,
+                          func: Callable,
+                          datasets: Dict[str, str],
+                          args: Dict[str, str] = {},
+                          constituents = None):
         """
         General function to add data on a regular grid.
 
@@ -1313,11 +1368,13 @@ class Result:
 
         at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings()
 
+        increments = self.place(list(datasets.values()),False)
+        if not increments: raise RuntimeError("received invalid dataset")
         with h5py.File(self.fname, 'a') as f:
-            for increment in self.place(datasets.values(),False).items():
+            for increment in increments.items():
                 for ty in increment[1].items():
                     for field in ty[1].items():
-                        d = list(field[1].values())[0]
+                        d: np.ma.MaskedArray = list(field[1].values())[0]
                         if np.any(d.mask): continue
                         dataset = {'f':{'data':np.reshape(d.data,tuple(self.cells)+d.data.shape[1:]),
                                         'label':list(datasets.values())[0],
@@ -1331,21 +1388,26 @@ class Result:
                                 result1 = result[at_cell_ho[x]]
 
                             path = '/'.join(['/',increment[0],ty[0],x,field[0]])
-                            dataset = f[path].create_dataset(r['label'],data=result1)
+                            h5_dataset = f[path].create_dataset(r['label'],data=result1)
 
                             now = datetime.datetime.now().astimezone()
-                            dataset.attrs['created'] = now.strftime('%Y-%m-%d %H:%M:%S%z') if h5py3 else \
-                                                       now.strftime('%Y-%m-%d %H:%M:%S%z').encode()
+                            h5_dataset.attrs['created'] = now.strftime('%Y-%m-%d %H:%M:%S%z') if h5py3 else \
+                                                          now.strftime('%Y-%m-%d %H:%M:%S%z').encode()
 
                             for l,v in r['meta'].items():
-                                dataset.attrs[l.lower()]=v if h5py3 else v.encode()
-                            creator = dataset.attrs['creator'] if h5py3 else \
-                                      dataset.attrs['creator'].decode()
-                            dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \
-                                                       f'damask.Result.{creator} v{damask.version}'.encode()
+                                h5_dataset.attrs[l.lower()]=v if h5py3 else v.encode()
+                            creator = h5_dataset.attrs['creator'] if h5py3 else \
+                                      h5_dataset.attrs['creator'].decode()
+                            h5_dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \
+                                                          f'damask.Result.{creator} v{damask.version}'.encode()
 
 
-    def _job_pointwise(self,group,func,datasets,args,lock):
+    def _job_pointwise(self,
+                       group: str,
+                       callback: Callable,
+                       datasets: Dict[str, str],
+                       args: Dict[str, str],
+                       lock: Lock) -> List[Union[None, Any]]:
         """Execute job for _add_generic_pointwise."""
         try:
             datasets_in = {}
@@ -1358,19 +1420,23 @@ class Result:
                                       'meta': {k:(v.decode() if not h5py3 and type(v) is bytes else v) \
                                                for k,v in loc.attrs.items()}}
             lock.release()
-            r = func(**datasets_in,**args)
+            r = callback(**datasets_in,**args)
             return [group,r]
         except Exception as err:
             print(f'Error during calculation: {err}.')
             return [None,None]
 
-    def _add_generic_pointwise(self,func,datasets,args={}):
+
+    def _add_generic_pointwise(self,
+                               func: Callable,
+                               datasets: Dict[str, Any],
+                               args: Dict[str, Any] = {}):
         """
         General function to add pointwise data.
 
         Parameters
         ----------
-        func : function
+        callback : function
             Callback function that calculates a new dataset from one or
             more datasets per HDF5 group.
         datasets : dictionary
@@ -1396,9 +1462,9 @@ class Result:
             print('No matching dataset found, no data was added.')
             return
 
-        default_arg = partial(self._job_pointwise,func=func,datasets=datasets,args=args,lock=lock)
+        default_arg = partial(self._job_pointwise,callback=func,datasets=datasets,args=args,lock=lock)
 
-        for group,result in util.show_progress(pool.imap_unordered(default_arg,groups),len(groups)):
+        for group,result in util.show_progress(pool.imap_unordered(default_arg,groups),len(groups)):# type: ignore
             if not result:
                 continue
             lock.acquire()
@@ -1410,15 +1476,14 @@ class Result:
                         dataset.attrs['overwritten'] = True
                     else:
                         shape = result['data'].shape
-                        if result['data'].size >= chunk_size*2:
+                        if compress := (result['data'].size >= chunk_size*2):
                             chunks = (chunk_size//np.prod(shape[1:]),)+shape[1:]
-                            compression = ('gzip',6)
                         else:
                             chunks = shape
-                            compression = (None,None)
                         dataset = f[group].create_dataset(result['label'],data=result['data'],
                                                           maxshape=shape, chunks=chunks,
-                                                          compression=compression[0], compression_opts=compression[1],
+                                                          compression = 'gzip' if compress else None,
+                                                          compression_opts = 6 if compress else None,
                                                           shuffle=True,fletcher32=True)
 
                     now = datetime.datetime.now().astimezone()
@@ -1440,7 +1505,8 @@ class Result:
         pool.join()
 
 
-    def export_XDMF(self,output='*'):
+    def export_XDMF(self,
+                    output: Union[str, List[str]] = '*'):
         """
         Write XDMF file to directly visualize data from DADF5 file.
 
@@ -1574,7 +1640,13 @@ class Result:
         return at_cell_ph,in_data_ph,at_cell_ho,in_data_ho
 
 
-    def export_VTK(self,output='*',mode='cell',constituents=None,fill_float=np.nan,fill_int=0,parallel=True):
+    def export_VTK(self,
+                   output: Union[str,list] = '*',
+                   mode: str = 'cell',
+                   constituents: IntSequence = None,
+                   fill_float: float = np.nan,
+                   fill_int: int = 0,
+                   parallel: bool = True):
         """
         Export to VTK cell/point data.
 
@@ -1612,12 +1684,12 @@ class Result:
         else:
             raise ValueError(f'invalid mode "{mode}"')
 
-        v.comments = util.execution_stamp('Result','export_VTK')
+        v.comments = [util.execution_stamp('Result','export_VTK')]
 
         N_digits = int(np.floor(np.log10(max(1,self.incs[-1]))))+1
 
         constituents_ = constituents if isinstance(constituents,Iterable) else \
-                      (range(self.N_constituents) if constituents is None else [constituents])
+                        (range(self.N_constituents) if constituents is None else [constituents])    # type: ignore
 
         suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \
                    [f'#{c}' for c in constituents_]
@@ -1637,7 +1709,7 @@ class Result:
 
                 for ty in ['phase','homogenization']:
                     for field in self.visible['fields']:
-                        outs = {}
+                        outs: Dict[str, np.ma.core.MaskedArray] = {}
                         for label in self.visible[ty+'s']:
                             if field not in f['/'.join([inc,ty,label])].keys(): continue
 
@@ -1665,7 +1737,10 @@ class Result:
                 v.save(f'{self.fname.stem}_inc{inc[10:].zfill(N_digits)}',parallel=parallel)
 
 
-    def get(self,output='*',flatten=True,prune=True):
+    def get(self,
+            output: Union[str, List[str]] = '*',
+            flatten: bool = True,
+            prune: bool = True):
         """
         Collect data per phase/homogenization reflecting the group/folder structure in the DADF5 file.
 
@@ -1687,7 +1762,7 @@ class Result:
             Datasets structured by phase/homogenization and according to selected view.
 
         """
-        r = {}
+        r = {}                                                                                      # type: ignore
 
         with h5py.File(self.fname,'r') as f:
             for inc in util.show_progress(self.visible['increments']):
@@ -1710,7 +1785,13 @@ class Result:
         return None if (type(r) == dict and r == {}) else r
 
 
-    def place(self,output='*',flatten=True,prune=True,constituents=None,fill_float=np.nan,fill_int=0):
+    def place(self,
+              output: Union[str, List[str]] = '*',
+              flatten: bool = True,
+              prune: bool = True,
+              constituents: IntSequence = None,
+              fill_float: float = np.nan,
+              fill_int: int = 0):
         """
         Merge data into spatial order that is compatible with the damask.VTK geometry representation.
 
@@ -1748,10 +1829,10 @@ class Result:
             Datasets structured by spatial position and according to selected view.
 
         """
-        r = {}
+        r = {}                                                                                      # type: ignore
 
-        constituents_ = constituents if isinstance(constituents,Iterable) else \
-                      (range(self.N_constituents) if constituents is None else [constituents])
+        constituents_ = list(map(int,constituents)) if isinstance(constituents,Iterable) else \
+                      (range(self.N_constituents) if constituents is None else [constituents])      # type: ignore
 
         suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \
                    [f'#{c}' for c in constituents_]
@@ -1797,7 +1878,9 @@ class Result:
         return None if (type(r) == dict and r == {}) else r
 
 
-    def export_setup(self,output='*',overwrite=False):
+    def export_setup(self,
+                     output: Union[str, List[str]] = '*',
+                     overwrite: bool = False):
         """
         Export configuration files.
 
@@ -1811,7 +1894,7 @@ class Result:
             Defaults to False.
 
         """
-        def export(name,obj,output,overwrite):
+        def export(name: str, obj: Union[h5py.Dataset,h5py.Group], output: Union[str,List[str]], overwrite: bool):
             if type(obj) == h5py.Dataset and  _match(output,[name]):
                 d = obj.attrs['description'] if h5py3 else obj.attrs['description'].decode()
                 if not Path(name).exists() or overwrite:
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 0e443ee62..61529eb46 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -994,7 +994,7 @@ class Rotation:
 
         """
         rng = np.random.default_rng(rng_seed)
-        r = rng.random(3 if shape is None else tuple(shape)+(3,) if hasattr(shape, '__iter__') else (shape,3)) #type: ignore
+        r = rng.random(3 if shape is None else tuple(shape)+(3,) if hasattr(shape, '__iter__') else (shape,3)) # type: ignore
 
         A = np.sqrt(r[...,2])
         B = np.sqrt(1.0-r[...,2])
diff --git a/python/damask/_table.py b/python/damask/_table.py
index 9ce0ce811..dc2c6f961 100644
--- a/python/damask/_table.py
+++ b/python/damask/_table.py
@@ -1,6 +1,6 @@
 import re
 import copy
-from typing import Union, Tuple, List
+from typing import Union, Tuple, List, Iterable
 
 import pandas as pd
 import numpy as np
@@ -12,9 +12,9 @@ class Table:
     """Manipulate multi-dimensional spreadsheet-like data."""
 
     def __init__(self,
-                 shapes: dict = {},
-                 data: np.ndarray = None,
-                 comments: Union[str, list] = None):
+                 shapes: dict,
+                 data: np.ndarray,
+                 comments: Union[str, Iterable[str]] = None):
         """
         New spreadsheet.
 
@@ -30,7 +30,7 @@ class Table:
 
         """
         comments_ = [comments] if isinstance(comments,str) else comments
-        self.comments = [] if comments_ is None else [c for c in comments_]
+        self.comments = [] if comments_ is None else [str(c) for c in comments_]
         self.shapes = { k:(v,) if isinstance(v,(np.int64,np.int32,int)) else v for k,v in shapes.items() }
         self.data = pd.DataFrame(data=data)
         self._relabel('uniform')
@@ -160,7 +160,7 @@ class Table:
             'linear'  ==> 1_v 2_v 3_v
 
         """
-        self.data.columns = self._label(self.shapes,how) #type: ignore
+        self.data.columns = self._label(self.shapes,how)                                            # type: ignore
 
 
     def _add_comment(self,
@@ -434,8 +434,8 @@ class Table:
 
 
     def rename(self,
-               old: Union[str, List[str]],
-               new: Union[str, List[str]],
+               old: Union[str, Iterable[str]],
+               new: Union[str, Iterable[str]],
                info: str = None) -> 'Table':
         """
         Rename column data.
diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index 6baa16214..c41470638 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -103,7 +103,7 @@ class VTK:
 
         Parameters
         ----------
-        comments : str or list of str
+        comments : str or sequence of str
             Comments.
 
         """
@@ -153,11 +153,11 @@ class VTK:
 
         Parameters
         ----------
-        cells : iterable of int, len (3)
+        cells : sequence of int, len (3)
             Number of cells along each dimension.
-        size : iterable of float, len (3)
+        size : sequence of float, len (3)
             Physical length along each dimension.
-        origin : iterable of float, len (3), optional
+        origin : sequence of float, len (3), optional
             Coordinates of grid origin.
 
         Returns
@@ -257,7 +257,7 @@ class VTK:
 
         Parameters
         ----------
-        grid : iterables of floats, len (3)
+        grid : sequence of sequences of floats, len (3)
             Grid coordinates along x, y, and z directions.
 
         Returns
diff --git a/python/tests/test_ConfigMaterial.py b/python/tests/test_ConfigMaterial.py
index d3ea986c5..f003e254e 100644
--- a/python/tests/test_ConfigMaterial.py
+++ b/python/tests/test_ConfigMaterial.py
@@ -113,15 +113,15 @@ class TestConfigMaterial:
     @pytest.mark.parametrize('N,n,kw',[
                                         (1,1,{'phase':'Gold',
                                               'O':[1,0,0,0],
-                                              'F_i':np.eye(3),
+                                              'V_e':np.eye(3),
                                               'homogenization':'SX'}),
                                         (3,1,{'phase':'Gold',
                                               'O':Rotation.from_random(3),
-                                              'F_i':np.broadcast_to(np.eye(3),(3,3,3)),
+                                              'V_e':np.broadcast_to(np.eye(3),(3,3,3)),
                                               'homogenization':'SX'}),
                                         (2,3,{'phase':np.broadcast_to(['a','b','c'],(2,3)),
                                               'O':Rotation.from_random((2,3)),
-                                              'F_i':np.broadcast_to(np.eye(3),(2,3,3,3)),
+                                              'V_e':np.broadcast_to(np.eye(3),(2,3,3,3)),
                                               'homogenization':['SX','PX']}),
                                         ])
     def test_material_add(self,kw,N,n):
diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90
index 86807507b..fcd2c189a 100644
--- a/src/HDF5_utilities.f90
+++ b/src/HDF5_utilities.f90
@@ -48,7 +48,7 @@ module HDF5_utilities
 !> @details for parallel IO, all dimension except for the last need to match
 !--------------------------------------------------------------------------------------------------
   interface HDF5_write
-#if defined(__GFORTRAN__) && __GNUC__<11
+#if defined(__GFORTRAN__)
     module procedure HDF5_write_real1
     module procedure HDF5_write_real2
     module procedure HDF5_write_real3
@@ -1214,7 +1214,7 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
 
 end subroutine HDF5_read_int7
 
-#if defined(__GFORTRAN__) && __GNUC__<11
+#if defined(__GFORTRAN__)
 
 !--------------------------------------------------------------------------------------------------
 !> @brief write dataset of type real with 1 dimension
@@ -1631,7 +1631,7 @@ subroutine HDF5_write_str(dataset,loc_id,datasetName)
 
 end subroutine HDF5_write_str
 
-#if defined(__GFORTRAN__) && __GNUC__<11
+#if defined(__GFORTRAN__)
 
 !--------------------------------------------------------------------------------------------------
 !> @brief write dataset of type integer with 1 dimension
diff --git a/src/IO.f90 b/src/IO.f90
index 240c9e992..047d11add 100644
--- a/src/IO.f90
+++ b/src/IO.f90
@@ -428,6 +428,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = 'too many systems requested'
     case (146)
       msg = 'number of values does not match'
+    case (147)
+      msg = 'V_e needs to be symmetric'
     case (148)
       msg = 'Nconstituents mismatch between homogenization and material'
 
diff --git a/src/material.f90 b/src/material.f90
index 2258c40a1..eafc411cd 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -27,7 +27,7 @@ module material
 
 
   type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0
-  type(tTensorContainer),   dimension(:), allocatable, public, protected :: material_F_i_0
+  type(tTensorContainer),   dimension(:), allocatable, public, protected :: material_V_e_0
 
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
@@ -133,7 +133,7 @@ subroutine parse()
   allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
 
   allocate(material_O_0(materials%length))
-  allocate(material_F_i_0(materials%length))
+  allocate(material_V_e_0(materials%length))
 
   allocate(ho_of(materials%length))
   allocate(ph_of(materials%length,homogenization_maxNconstituents),source=-1)
@@ -154,7 +154,7 @@ subroutine parse()
         if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
 
         allocate(material_O_0(ma)%data(constituents%length))
-        allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length))
+        allocate(material_V_e_0(ma)%data(1:3,1:3,constituents%length))
 
         do co = 1, constituents%length
           constituent => constituents%get(co)
@@ -162,7 +162,9 @@ subroutine parse()
           ph_of(ma,co) = phases%getIndex(constituent%get_asString('phase'))
 
           call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
-          material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3])
+          material_V_e_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('V_e',defaultVal=math_I3,requiredShape=[3,3])
+          if (any(dNeq(material_V_e_0(ma)%data(1:3,1:3,co),transpose(material_V_e_0(ma)%data(1:3,1:3,co))))) &
+            call IO_error(147)
 
         end do
         if (dNeq(sum(v_of(ma,:)),1.0_pReal,1.e-9_pReal)) call IO_error(153,ext_msg='constituent')
diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index fc9e40e02..88b86a8d9 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -210,7 +210,6 @@ module subroutine mechanical_init(phases)
     Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
-
     phase, &
     mech
 
@@ -239,11 +238,8 @@ module subroutine mechanical_init(phases)
 
     allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
@@ -265,23 +261,19 @@ module subroutine mechanical_init(phases)
     ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
     do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
       ph = material_phaseID(co,ce)
-      phase_mechanical_Fi0(ph)%data(1:3,1:3,material_phaseEntry(co,ce)) = material_F_i_0(ma)%data(1:3,1:3,co)
+      en = material_phaseEntry(co,ce)
+      phase_mechanical_F(ph)%data(1:3,1:3,en)  = math_I3
+      phase_mechanical_Fp(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%asMatrix()              ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
+      phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), &
+                                                        transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
+      phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co)))
     enddo
   enddo
 
   do ph = 1, phases%length
-    do en = 1, count(material_phaseID == ph)
-
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_O_0(ph)%data(en)%asMatrix()                 ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) &
-                                                / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal)
-      phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3
-      phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), &
-                                                                   phase_mechanical_Fp0(ph)%data(1:3,1:3,en)))  ! assuming that euler angles are given in internal strain free configuration
-    enddo
-    phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data
-    phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data
-    phase_mechanical_F(ph)%data  = phase_mechanical_F0(ph)%data
+    phase_mechanical_F0(ph)%data  = phase_mechanical_F(ph)%data
+    phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data
+    phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data
   enddo
 
 
@@ -318,7 +310,6 @@ module subroutine mechanical_init(phases)
 
   end select
 
-
   call eigen_init(phases)
 
 
@@ -1244,6 +1235,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
 end function phase_mechanical_dPdF
 
 
+!--------------------------------------------------------------------------------------------------
+!< @brief Write restart information to file.
+!--------------------------------------------------------------------------------------------------
 module subroutine mechanical_restartWrite(groupHandle,ph)
 
   integer(HID_T), intent(in) :: groupHandle
@@ -1251,36 +1245,41 @@ module subroutine mechanical_restartWrite(groupHandle,ph)
 
 
   call HDF5_write(plasticState(ph)%state,groupHandle,'omega_plastic')
-  call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
-  call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
-  call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
-  call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p')
   call HDF5_write(phase_mechanical_S(ph)%data,groupHandle,'S')
   call HDF5_write(phase_mechanical_F(ph)%data,groupHandle,'F')
+  call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p')
+  call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
+  call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
+  call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
 
 end subroutine mechanical_restartWrite
 
 
+!--------------------------------------------------------------------------------------------------
+!< @brief Read restart information from file.
+!--------------------------------------------------------------------------------------------------
 module subroutine mechanical_restartRead(groupHandle,ph)
 
   integer(HID_T), intent(in) :: groupHandle
   integer, intent(in) :: ph
 
+  integer :: en
+
 
   call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic')
-  call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
-  call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
-  call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')
-  call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p')
   call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S')
   call HDF5_read(phase_mechanical_F0(ph)%data,groupHandle,'F')
+  call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p')
+  call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
+  call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')
+  call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
 
 end subroutine mechanical_restartRead
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
+!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics).
+!--------------------------------------------------------------------------------------------------
 module function mechanical_S(ph,en) result(S)
 
   integer, intent(in) :: ph,en
@@ -1292,9 +1291,9 @@ module function mechanical_S(ph,en) result(S)
 end function mechanical_S
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Get plastic velocity gradient (for use by non-mech physics)
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
+!< @brief Get plastic velocity gradient (for use by non-mech physics).
+!--------------------------------------------------------------------------------------------------
 module function mechanical_L_p(ph,en) result(L_p)
 
   integer, intent(in) :: ph,en
@@ -1306,9 +1305,9 @@ module function mechanical_L_p(ph,en) result(L_p)
 end function mechanical_L_p
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Get elastic deformation gradient (for use by non-mech physics)
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
+!< @brief Get elastic deformation gradient (for use by non-mech physics).
+!--------------------------------------------------------------------------------------------------
 module function mechanical_F_e(ph,en) result(F_e)
 
   integer, intent(in) :: ph,en
@@ -1320,9 +1319,9 @@ module function mechanical_F_e(ph,en) result(F_e)
 end function mechanical_F_e
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Get second Piola-Kichhoff stress (for use by homogenization)
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
+!< @brief Get second Piola-Kichhoff stress (for use by homogenization).
+!--------------------------------------------------------------------------------------------------
 module function phase_P(co,ce) result(P)
 
   integer, intent(in) :: co, ce
@@ -1334,9 +1333,9 @@ module function phase_P(co,ce) result(P)
 end function phase_P
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Get deformation gradient (for use by homogenization)
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
+!< @brief Get deformation gradient (for use by homogenization).
+!--------------------------------------------------------------------------------------------------
 module function phase_F(co,ce) result(F)
 
   integer, intent(in) :: co, ce
@@ -1348,9 +1347,9 @@ module function phase_F(co,ce) result(F)
 end function phase_F
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Set deformation gradient (for use by homogenization)
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
+!< @brief Set deformation gradient (for use by homogenization).
+!--------------------------------------------------------------------------------------------------
 module subroutine phase_set_F(F,co,ce)
 
   real(pReal), dimension(3,3), intent(in) :: F
diff --git a/src/rotations.f90 b/src/rotations.f90
index 7aa81fcab..5367f2c07 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -323,17 +323,16 @@ pure function rotTensor2(self,T,active) result(tRot)
 
   logical           :: passive
 
+
   if (present(active)) then
     passive = .not. active
   else
     passive = .true.
   endif
 
-  if (passive) then
-    tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix()))
-  else
-    tRot = matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix())
-  endif
+  tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
+               matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
+               passive)
 
 end function rotTensor2
 
@@ -450,6 +449,7 @@ pure function qu2om(qu) result(om)
   om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
 
   if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
+  om = om/math_det33(om)**(1.0_pReal/3.0_pReal)
 
 end function qu2om
 
@@ -575,7 +575,7 @@ end function qu2cu
 
 !---------------------------------------------------------------------------------------------------
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief convert rotation matrix to cubochoric
+!> @brief convert rotation matrix to unit quaternion
 !> @details the original formulation (direct conversion) had (numerical?) issues
 !---------------------------------------------------------------------------------------------------
 pure function om2qu(om) result(qu)
@@ -602,14 +602,15 @@ pure function om2qu(om) result(qu)
       endif
   endif
   if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
-  qu = qu*[1.0_pReal,P,P,P]
+  qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
+  qu = qu/norm2(qu)
 
 end function om2qu
 
 
 !---------------------------------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
-!> @brief orientation matrix to Euler angles
+!> @brief convert orientation matrix to Euler angles
 !> @details Two step check for special cases to avoid invalid operations (not needed for python)
 !---------------------------------------------------------------------------------------------------
 pure function om2eu(om) result(eu)
@@ -1334,8 +1335,8 @@ pure function cu2ho(cu) result(ho)
       ! transform to sphere grid (inverse Lambert)
       ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
       c = sum(T**2)
-      s = PI * c/(24.0*XYZ(3)**2)
-      c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3)
+      s = c * PI/(24.0*XYZ(3)**2)
+      c = c * sqrt(PI/24.0_pReal) / XYZ(3)
       q = sqrt( 1.0 - s )
       LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
     end if special