diff --git a/VERSION b/VERSION
index 2b4efa3f2..9bb6a55a4 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-3.0.0-alpha7-753-gaa28fe677
+3.0.0-alpha7-784-g69239661c
diff --git a/python/damask/_crystal.py b/python/damask/_crystal.py
index 9434515f5..adcaf610e 100644
--- a/python/damask/_crystal.py
+++ b/python/damask/_crystal.py
@@ -2,11 +2,11 @@ from typing import Optional, Union, Dict, List, Tuple
 
 import numpy as np
 
-from ._typehints import FloatSequence, CrystalFamily, CrystalLattice, CrystalKinematics
+from ._typehints import FloatSequence, CrystalFamily, BravaisLattice, CrystalKinematics
 from . import util
 from . import Rotation
 
-lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
+lattice_symmetries: Dict[BravaisLattice, CrystalFamily] = {
                 'aP': 'triclinic',
 
                 'mP': 'monoclinic',
@@ -27,7 +27,7 @@ lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
                 'cF': 'cubic',
                }
 
-orientation_relationships: Dict[str, Dict[CrystalLattice,np.ndarray]] = {
+orientation_relationships: Dict[str, Dict[BravaisLattice,np.ndarray]] = {
  'KS': {
    'cF': np.array([
        [[-1, 0, 1],[ 1, 1, 1]],
@@ -323,7 +323,7 @@ class Crystal():
 
     def __init__(self, *,
                  family: Optional[CrystalFamily] = None,
-                 lattice: Optional[CrystalLattice] = None,
+                 lattice: Optional[BravaisLattice] = None,
                  a: Optional[float] = None, b: Optional[float] = None, c: Optional[float] = None,
                  alpha: Optional[float] = None, beta: Optional[float] = None, gamma: Optional[float] = None,
                  degrees: bool = False):
@@ -548,7 +548,17 @@ class Crystal():
 
     @property
     def symmetry_operations(self) -> Rotation:
-        """Symmetry operations as Rotations."""
+        """
+        Return symmetry operations.
+
+        References
+        ----------
+        U.F. Kocks et al.,
+        Texture and Anisotropy:
+        Preferred Orientations in Polycrystals and their Effect on Materials Properties.
+        Cambridge University Press 1998. Table II
+
+        """
         _symmetry_operations: Dict[CrystalFamily, List]  = {
             'cubic':         [
                               [ 1.0,            0.0,            0.0,            0.0            ],
@@ -772,7 +782,7 @@ class Crystal():
             Directions and planes of deformation mode families.
 
         """
-        _kinematics: Dict[CrystalLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
+        _kinematics: Dict[BravaisLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
             'cF': {
                 'slip': [np.array([
                            [+0,+1,-1, +1,+1,+1],
@@ -1015,7 +1025,7 @@ class Crystal():
 
 
     def relation_operations(self,
-                            model: str) -> Tuple[CrystalLattice, Rotation]:
+                            model: str) -> Tuple[BravaisLattice, Rotation]:
         """
         Crystallographic orientation relationships for phase transformations.
 
diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index 4db24bae5..cf6684458 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -3,7 +3,7 @@ from typing import Optional, Union, TypeVar
 
 import numpy as np
 
-from ._typehints import FloatSequence, IntSequence, CrystalFamily, CrystalLattice
+from ._typehints import FloatSequence, IntSequence, CrystalFamily, BravaisLattice
 from . import Rotation
 from . import Crystal
 from . import util
@@ -73,7 +73,7 @@ class Orientation(Rotation,Crystal):
                  rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
                  *,
                  family: Optional[CrystalFamily] = None,
-                 lattice: Optional[CrystalLattice] = None,
+                 lattice: Optional[BravaisLattice] = None,
                  a: Optional[float] = None, b: Optional[float] = None, c: Optional[float] = None,
                  alpha: Optional[float] = None, beta: Optional[float] = None, gamma: Optional[float] = None,
                  degrees: bool = False):
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 118234bf1..cd3164db4 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -307,7 +307,8 @@ class Rotation:
             p_m = self.quaternion[...,1:]
             q_o = other.quaternion[...,0:1]
             p_o = other.quaternion[...,1:]
-            q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,)))
+            qmo = q_m*q_o
+            q = (qmo - np.einsum('...i,...i',p_m,p_o).reshape(qmo.shape))
             p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
             return self.copy(Rotation(np.block([q,p]))._standardize())
         else:
@@ -1323,28 +1324,41 @@ class Rotation:
 
     @staticmethod
     def _qu2eu(qu: np.ndarray) -> np.ndarray:
-        """Quaternion to Bunge Euler angles."""
-        q02   = qu[...,0:1]*qu[...,2:3]
-        q13   = qu[...,1:2]*qu[...,3:4]
-        q01   = qu[...,0:1]*qu[...,1:2]
-        q23   = qu[...,2:3]*qu[...,3:4]
-        q03_s = qu[...,0:1]**2+qu[...,3:4]**2
-        q12_s = qu[...,1:2]**2+qu[...,2:3]**2
-        chi = np.sqrt(q03_s*q12_s)
+        """
+        Quaternion to Bunge Euler angles.
 
-        eu = np.where(np.abs(q12_s) < 1.e-8,
-                np.block([np.arctan2(-_P*2.*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2),
-                          np.zeros(qu.shape[:-1]+(2,))]),
-                      np.where(np.abs(q03_s) < 1.e-8,
-                          np.block([np.arctan2(   2.*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
-                                    np.broadcast_to(np.pi,qu[...,0:1].shape),
-                                    np.zeros(qu.shape[:-1]+(1,))]),
-                          np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi),
-                                    np.arctan2(           2.*chi,    q03_s-q12_s   ),
-                                    np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)])
-                              )
-                     )
-        eu[np.abs(eu) < 1.e-6] = 0.
+        References
+        ----------
+        E. Bernardes and S. Viollet, PLoS ONE 17(11):e0276302, 2022
+        https://doi.org/10.1371/journal.pone.0276302
+
+        """
+        a =     qu[...,0:1]
+        b = -_P*qu[...,3:4]
+        c = -_P*qu[...,1:2]
+        d = -_P*qu[...,2:3]
+
+        eu = np.block([
+            np.arctan2(b,a),
+            np.arccos(2*(a**2+b**2)/(a**2+b**2+c**2+d**2)-1),
+            np.arctan2(-d,c),
+        ])
+
+        eu_sum  = eu[...,0] + eu[...,2]
+        eu_diff = eu[...,0] - eu[...,2]
+
+        is_zero  = np.isclose(eu[...,1],0.0)
+        is_pi    = np.isclose(eu[...,1],np.pi)
+        is_ok    = ~np.logical_or(is_zero,is_pi)
+
+        eu[...,0][is_zero] =  2*eu[...,0][is_zero]
+        eu[...,0][is_pi]   = -2*eu[...,2][is_pi]
+        eu[...,2][~is_ok]  = 0.0
+        eu[...,0][is_ok]   = eu_diff[is_ok]
+        eu[...,2][is_ok]   = eu_sum [is_ok]
+
+        eu[np.logical_or(np.abs(eu)         < 1.e-6,
+                         np.abs(eu-2*np.pi) < 1.e-6)] = 0.
         return np.where(eu < 0., eu%(np.pi*np.array([2.,1.,2.])),eu)
 
     @staticmethod
diff --git a/python/damask/_typehints.py b/python/damask/_typehints.py
index 5891e6305..29ab9f3ad 100644
--- a/python/damask/_typehints.py
+++ b/python/damask/_typehints.py
@@ -11,7 +11,7 @@ IntSequence = Union[np.ndarray,Sequence[int]]
 StrSequence = Union[np.ndarray,Sequence[str]]
 FileHandle = Union[TextIO, str, Path]
 CrystalFamily = Union[None,Literal['triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic']]
-CrystalLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
+BravaisLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
 CrystalKinematics = Literal['slip', 'twin']
 NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
 # BitGenerator does not exists in older numpy versions
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 6e135212f..6b9cd2b1a 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -975,6 +975,13 @@ class TestRotation:
             assert np.allclose(rot_broadcast.quaternion[...,i,:], rot.quaternion)
 
 
+    @pytest.mark.parametrize('shape',[(3,2),(4,6)])
+    def test_broadcastcomposition(self,shape):
+        a = Rotation.from_random(shape[0])
+        b = Rotation.from_random(shape[1])
+        assert (a[:,np.newaxis]*b[np.newaxis,:]).allclose(a.broadcast_to(shape)*b.broadcast_to(shape))
+
+
     @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion,        np.array([-1,0,0,0])),
                                                  (Rotation.from_quaternion,        np.array([1,1,1,0])),
                                                  (Rotation.from_Euler_angles,      np.array([1,4,0])),
diff --git a/src/crystal.f90 b/src/crystal.f90
index ae460a811..7bf9dc2f7 100644
--- a/src/crystal.f90
+++ b/src/crystal.f90
@@ -1902,7 +1902,7 @@ end function buildInteraction
 !> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
 !> @details Order: Direction, plane (normal), and common perpendicular
 !--------------------------------------------------------------------------------------------------
-function buildCoordinateSystem(active,potential,system,lattice,cOverA)
+function buildCoordinateSystem(active,potential,system,lattice,cOverA) result(coordinateSystem)
 
   integer, dimension(:), intent(in) :: &
     active, &                                                                                       !< # of active systems per family
@@ -1914,7 +1914,7 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
   real(pREAL),                 intent(in) :: &
     cOverA
   real(pREAL), dimension(3,3,sum(active)) :: &
-    buildCoordinateSystem
+    coordinateSystem
 
   real(pREAL), dimension(3) :: &
     direction, normal
@@ -1937,10 +1937,14 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
 
       select case(lattice)
 
-        case ('cF','cI','tI')
+        case ('cF','cI')
           direction = system(1:3,p)
           normal    = system(4:6,p)
 
+        case ('tI')
+          direction = [ system(1,p), system(2,p), system(3,p)*cOverA ]
+          normal    = [ system(4,p), system(5,p), system(6,p)/cOverA ]
+
         case ('hP')
           direction = [ system(1,p)*1.5_pREAL, &
                        (system(1,p)+2.0_pREAL*system(2,p))*sqrt(0.75_pREAL), &
@@ -1954,10 +1958,10 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
 
       end select
 
-      buildCoordinateSystem(1:3,1,a) = direction/norm2(direction)
-      buildCoordinateSystem(1:3,2,a) = normal   /norm2(normal)
-      buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
-                                                  normal   /norm2(normal))
+      coordinateSystem(1:3,1,a) = direction/norm2(direction)
+      coordinateSystem(1:3,2,a) = normal   /norm2(normal)
+      coordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
+                                             normal   /norm2(normal))
 
     end do activeSystems
   end do activeFamilies
@@ -2245,6 +2249,12 @@ subroutine crystal_selfTest
   CoSy   = buildCoordinateSystem([1],[1],system,'cF',0.0_pREAL)
   if (any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem'
 
+  if (any(dNeq(buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'cI',0.0_pReal), &
+               buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'tI',1.0_pReal)))) &
+                                                           error stop 'cI/tI coordinate system'
+  if (all(dEq( buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'tI',1.1_pReal + r(1)*0.9_pReal), &
+               buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'tI',1.0_pReal)))) &
+                                                           error stop 'tI coordinate system'
   do i = 1, 10
     call random_number(C)
     C_cF = crystal_symmetrize_C66(C,'cI')
diff --git a/src/phase.f90 b/src/phase.f90
index 4b01c6ef0..4cc341913 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -427,12 +427,12 @@ subroutine phase_init
     phase => phases%get_dict(ph)
     refs = config_listReferences(phase,indent=3)
     if (len(refs) > 0) print'(/,1x,a)', refs
+    phase_rho(ph) = phase%get_asReal('rho',defaultVal=0.0_pREAL)
     phase_lattice(ph) = phase%get_asStr('lattice')
     if (all(phase_lattice(ph) /= ['cF','cI','hP','tI'])) &
       call IO_error(130,ext_msg='phase_init: '//phase%get_asStr('lattice'))
     if (any(phase_lattice(ph) == ['hP','tI'])) &
       phase_cOverA(ph) = phase%get_asReal('c/a')
-    phase_rho(ph) = phase%get_asReal('rho',defaultVal=0.0_pREAL)
     allocate(phase_O_0(ph)%data(count(material_ID_phase==ph)))
   end do
 
diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90
index e2ec9922c..fc3fbf5a2 100644
--- a/src/phase_damage_anisobrittle.f90
+++ b/src/phase_damage_anisobrittle.f90
@@ -50,7 +50,7 @@ module function anisobrittle_init() result(mySources)
   if (count(mySources) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:damage:anisobrittle init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
@@ -64,7 +64,7 @@ module function anisobrittle_init() result(mySources)
 
       associate(prm  => param(ph))
 
-        print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
+        print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
         refs = config_listReferences(src,indent=3)
         if (len(refs) > 0) print'(/,1x,a)', refs
 
diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90
index 2efff9f7d..39bb762cf 100644
--- a/src/phase_damage_isobrittle.f90
+++ b/src/phase_damage_isobrittle.f90
@@ -48,7 +48,7 @@ module function isobrittle_init() result(mySources)
   if (count(mySources) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:damage:isobrittle init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
@@ -66,7 +66,7 @@ module function isobrittle_init() result(mySources)
 
         prm%W_crit = src%get_asReal('G_crit')/src%get_asReal('l_c')
 
-        print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
+        print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
         refs = config_listReferences(src,indent=3)
         if (len(refs) > 0) print'(/,1x,a)', refs
 
diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90
index 5bdb57f57..8e0f7209a 100644
--- a/src/phase_mechanical_eigen.f90
+++ b/src/phase_mechanical_eigen.f90
@@ -95,34 +95,6 @@ function kinematics_active(kinematics_label,kinematics_length)  result(active_ki
 end function kinematics_active
 
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Checks if a damage kinematic mechanism is active.
-!--------------------------------------------------------------------------------------------------
-function kinematics_active2(kinematics_label)  result(active_kinematics)
-
-  character(len=*), intent(in)       :: kinematics_label                                            !< name of kinematic mechanism
-  logical, dimension(:), allocatable :: active_kinematics
-
-  type(tDict), pointer :: &
-    phases, &
-    phase, &
-    kinematics_type
-  integer :: ph
-
-
-  phases => config_material%get_dict('phase')
-  allocate(active_kinematics(phases%length), source = .false.)
-  do ph = 1, phases%length
-    phase => phases%get_dict(ph)
-    kinematics_type => phase%get_dict('damage',defaultVal=emptyDict)
-    active_kinematics(ph) = kinematics_type%get_asStr('type',defaultVal='n/a') == kinematics_label
-  end do
-
-
-end function kinematics_active2
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient
 ! ToDo: MD: S is Mi?
diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90
index 027f71c81..5bc64bb49 100644
--- a/src/phase_mechanical_eigen_thermalexpansion.f90
+++ b/src/phase_mechanical_eigen_thermalexpansion.f90
@@ -27,7 +27,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
   integer, intent(in)                  :: kinematics_length
   logical, dimension(:,:), allocatable :: myKinematics
 
-  integer :: Ninstances, p, k
+  integer :: p, k
   type(tList), pointer :: &
     kinematics
   type(tDict), pointer :: &
@@ -37,15 +37,13 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
 
 
   myKinematics = kinematics_active('thermalexpansion',kinematics_length)
-  Ninstances = count(myKinematics)
-  print'(/,a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
-  if (Ninstances == 0) return
+  if (count(myKinematics) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:eigen:thermalexpansion init  -+>>>'
-
+  print'(/,1x,a,1x,i0)', '# phases:',count(myKinematics); flush(IO_STDOUT)
 
   phases => config_material%get_dict('phase')
-  allocate(param(Ninstances))
+  allocate(param(count(myKinematics)))
   allocate(kinematics_thermal_expansion_instance(phases%length), source=0)
 
   do p = 1, phases%length
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index 89b5027ac..40b5e4cf1 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -34,7 +34,7 @@ module subroutine elastic_init(phases)
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic init  -+>>>'
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic:Hooke init  -+>>>'
 
-  print'(/,a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',phases%length; flush(IO_STDOUT)
 
 
   allocate(param(phases%length))
@@ -43,7 +43,7 @@ module subroutine elastic_init(phases)
     phase   => phases%get_dict(ph)
     mech    => phase%get_dict('mechanical')
     elastic => mech%get_dict('elastic')
-    print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
+    print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
     refs = config_listReferences(elastic,indent=3)
     if (len(refs) > 0) print'(/,1x,a)', refs
     if (elastic%get_asStr('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asStr('type'))
diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90
index 08dd3decd..a3567efe4 100644
--- a/src/phase_mechanical_plastic_dislotungsten.f90
+++ b/src/phase_mechanical_plastic_dislotungsten.f90
@@ -108,11 +108,11 @@ module function plastic_dislotungsten_init() result(myPlasticity)
   if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotungsten init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
   print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
 
+  print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get_dict('phase')
   allocate(param(phases%length))
diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index a1f22baa7..a93c4752a 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -155,7 +155,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
   if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
   print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'
@@ -166,6 +165,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
   print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
   print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032'
 
+  print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get_dict('phase')
   allocate(param(phases%length))
diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90
index be859dddb..e7935ddf1 100644
--- a/src/phase_mechanical_plastic_isotropic.f90
+++ b/src/phase_mechanical_plastic_isotropic.f90
@@ -68,11 +68,11 @@ module function plastic_isotropic_init() result(myPlasticity)
   if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   print'(/,1x,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
   print'(  1x,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
 
+  print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get_dict('phase')
   allocate(param(phases%length))
diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90
index 2ce34cd95..81fa188c8 100644
--- a/src/phase_mechanical_plastic_kinehardening.f90
+++ b/src/phase_mechanical_plastic_kinehardening.f90
@@ -100,11 +100,11 @@ module function plastic_kinehardening_init() result(myPlasticity)
   if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   print'(/,1x,a)', 'J.A. Wollmershauser et al., International Journal of Fatigue 36:181–193, 2012'
   print'(  1x,a)', 'https://doi.org/10.1016/j.ijfatigue.2011.07.008'
 
+  print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get_dict('phase')
   allocate(param(phases%length))
diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90
index 748e0579b..b42c0513c 100644
--- a/src/phase_mechanical_plastic_none.f90
+++ b/src/phase_mechanical_plastic_none.f90
@@ -25,7 +25,7 @@ module function plastic_none_init() result(myPlasticity)
   if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:none init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90
index 3e34256f1..e90dcae3c 100644
--- a/src/phase_mechanical_plastic_nonlocal.f90
+++ b/src/phase_mechanical_plastic_nonlocal.f90
@@ -203,7 +203,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
   end if
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:nonlocal init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
 
   print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
   print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'
@@ -211,6 +210,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
   print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
   print'(  1x,a)', 'http://publications.rwth-aachen.de/record/229993'
 
+  print'(/,1x,a,1x,i0)', '# phases:',Ninstances; flush(IO_STDOUT)
 
   phases => config_material%get_dict('phase')
   allocate(geom(phases%length))
diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90
index 8a8018892..4091e0f76 100644
--- a/src/phase_mechanical_plastic_phenopowerlaw.f90
+++ b/src/phase_mechanical_plastic_phenopowerlaw.f90
@@ -105,7 +105,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:phenopowerlaw init  -+>>>'
-  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
@@ -124,7 +124,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     mech => phase%get_dict('mechanical')
     pl => mech%get_dict('plastic')
 
-    print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
+    print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
     refs = config_listReferences(pl,indent=3)
     if (len(refs) > 0) print'(/,1x,a)', refs
 
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index e21f173ad..f63eb5e19 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -76,7 +76,9 @@ module subroutine thermal_init(phases)
     thermal
   type(tList), pointer :: &
     sources
-  character(len=:), allocatable :: refs
+  character(len=:), allocatable :: &
+    refs, &
+    extmsg
   integer :: &
     ph, so, &
     Nmembers
@@ -88,6 +90,7 @@ module subroutine thermal_init(phases)
   allocate(thermalState(phases%length))
   allocate(thermal_Nsources(phases%length),source = 0)
   allocate(param(phases%length))
+  extmsg = ''
 
   do ph = 1, phases%length
     Nmembers = count(material_ID_phase == ph)
@@ -106,6 +109,15 @@ module subroutine thermal_init(phases)
       if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33')
       param(ph)%K = crystal_symmetrize_33(param(ph)%K,phase_lattice(ph))
 
+      ! sanity checks
+      if (    param(ph)%C_p <= 0.0_pREAL )  extmsg = trim(extmsg)//' C_p'
+      if (any(param(ph)%K   <  0.0_pREAL))  extmsg = trim(extmsg)//' K'
+      if (    phase_rho(ph) <= 0.0_pREAL )  extmsg = trim(extmsg)//' rho'
+
+!--------------------------------------------------------------------------------------------------
+!  exit if any parameter is out of range
+      if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
+
 #if defined(__GFORTRAN__)
       param(ph)%output = output_as1dStr(thermal)
 #else
diff --git a/src/phase_thermal_source_dissipation.f90 b/src/phase_thermal_source_dissipation.f90
index 760259286..f0e75b10a 100644
--- a/src/phase_thermal_source_dissipation.f90
+++ b/src/phase_thermal_source_dissipation.f90
@@ -42,7 +42,7 @@ module function source_dissipation_init(maxNsources) result(isMySource)
   if (count(isMySource) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:thermal:source_dissipation init  -+>>>'
-  print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',count(isMySource); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
@@ -60,7 +60,7 @@ module function source_dissipation_init(maxNsources) result(isMySource)
       if (isMySource(so,ph)) then
         associate(prm  => param(ph))
           src => sources%get_dict(so)
-          print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
+          print'(/,1x,a,1x,i0,1x,a,1x,a,1x,i0)', 'phase',ph,'('//phases%key(ph)//')','source',so
           refs = config_listReferences(src,indent=3)
           if (len(refs) > 0) print'(/,1x,a)', refs
 
diff --git a/src/phase_thermal_source_externalheat.f90 b/src/phase_thermal_source_externalheat.f90
index 1a803549c..bd642af48 100644
--- a/src/phase_thermal_source_externalheat.f90
+++ b/src/phase_thermal_source_externalheat.f90
@@ -44,7 +44,7 @@ module function source_externalheat_init(maxNsources) result(isMySource)
   if (count(isMySource) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:thermal:source_externalheat init  -+>>>'
-  print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
+  print'(/,1x,a,1x,i0)', '# phases:',count(isMySource); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
@@ -64,7 +64,7 @@ module function source_externalheat_init(maxNsources) result(isMySource)
         source_ID(ph) = so
         associate(prm  => param(ph))
           src => sources%get_dict(so)
-          print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
+          print'(/,1x,a,1x,i0,1x,a,1x,a,1x,i0)', 'phase',ph,'('//phases%key(ph)//')','source',so
           refs = config_listReferences(src,indent=3)
           if (len(refs) > 0) print'(/,1x,a)', refs