Merge branch 'misc-improvements' into 'development'

Misc improvements See merge request damask/DAMASK!54
2019-01-05 15:20:07 +01:00 · 2019-01-05 15:20:07 +01:00 · 53bc24cced
parent 3c5df0a4a4 45d11f81b0
commit 53bc24cced
56 changed files with 544 additions and 4424 deletions
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@ -158,12 +158,12 @@ Post_AverageDown:
    - master
    - release

-#Post_General:
-#  stage: postprocessing
-#  script: PostProcessing/test.py
-#  except:
-#    - master
-#    - release
+Post_General:
+  stage: postprocessing
+  script: PostProcessing/test.py
+  except:
+    - master
+    - release

 Post_GeometryReconstruction:
  stage: postprocessing
@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip:
    - master
    - release

-#TextureComponents:
-#  stage: spectral
-#  script: TextureComponents/test.py
-#  except:
-#    - master
-#    - release
+TextureComponents:
+  stage: spectral
+  script: TextureComponents/test.py
+  except:
+    - master
+    - release


 ###################################################################################################
@ -468,27 +468,24 @@ AbaqusStd:
  script:
    - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
    - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT abaqus
-  except:
-    - master
-    - release
+  only:
+    - development

 Marc:
  stage: createDocumentation
  script:
    - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
    - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT marc
-  except:
-    - master
-    - release
+  only:
+    - development

 Spectral:
  stage: createDocumentation
  script:
    - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
    - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT spectral
-  except:
-    - master
-    - release
+  only:
+    - development

 ##################################################################################################
 backupData:
--- a/4
+++ b/4
@ -7,11 +7,11 @@ all: spectral FEM processing

 .PHONY: spectral
 spectral: build/spectral
-	@(cd build/spectral;make --no-print-directory -ws all install;)
+	@(cd build/spectral;make -j4 --no-print-directory -ws all install;)

 .PHONY: FEM
 FEM: build/FEM
-	@(cd build/FEM; make --no-print-directory -ws all install;)
+	@(cd build/FEM; make -j4 --no-print-directory -ws all install;)

 .PHONY: build/spectral
 build/spectral:
--- a/2
+++ b/2
@ -1 +1 @@
-Subproject commit b9a52a85cd65cc27a8e863302bd984abdcad1455
+Subproject commit 59b0cbe899f272476fb6f00f0f8860428e6ceba3
--- a/env/DAMASK.csh
+++ b/env/DAMASK.csh
@ -64,7 +64,7 @@ endif

 setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS
 if ( ! $?PYTHONPATH ) then
-  setenv PYTHONPATH $DAMASK_ROOT/lib
+  setenv PYTHONPATH $DAMASK_ROOT/python
 else
-  setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH
+  setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH
 endif
--- a/env/DAMASK.sh
+++ b/env/DAMASK.sh
@ -95,7 +95,7 @@ if [ ! -z "$PS1" ]; then
 fi

 export DAMASK_NUM_THREADS
-export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH
+export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH

 for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do
  unset "${var}"
--- a/env/DAMASK.zsh
+++ b/env/DAMASK.zsh
@ -88,7 +88,7 @@ if [ ! -z "$PS1" ]; then
 fi

 export DAMASK_NUM_THREADS
-export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH
+export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH

 for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do
  unset "${var}"
--- a/processing/post/addEuclideanDistance.py
+++ b/processing/post/addEuclideanDistance.py
@ -149,7 +149,6 @@ for name in filenames:

  errors  = []
  remarks = []
-  column = {}
  
  if not 3 >= table.label_dimension(options.pos) >= 1:
    errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
--- a/processing/post/addMises.py
+++ b/processing/post/addMises.py
@ -4,6 +4,7 @@
 import os,sys,math
 import numpy as np
 from optparse import OptionParser
+from collections import OrderedDict
 import damask

 scriptName = os.path.splitext(os.path.basename(__file__))[0]
@ -63,10 +64,10 @@ for name in filenames:

 # ------------------------------------------ sanity checks ----------------------------------------

-  items = {
-            'strain': {'dim': 9, 'shape': [3,3], 'labels':options.strain, 'active':[], 'column': []},
-            'stress': {'dim': 9, 'shape': [3,3], 'labels':options.stress, 'active':[], 'column': []},
-          }
+  items = OrderedDict([
+            ('strain', {'dim': 9, 'shape': [3,3], 'labels':options.strain, 'active':[], 'column': []}),
+            ('stress', {'dim': 9, 'shape': [3,3], 'labels':options.stress, 'active':[], 'column': []})
+          ])
  errors  = []
  remarks = []
  
--- a/processing/post/postResults.py
+++ b/processing/post/postResults.py
@ -830,10 +830,9 @@ if options.info:

 elementsOfNode = {}
 Nelems = stat['NumberOfElements']
-starttime = time.time()
 for e in range(Nelems):
-  if options.verbose and Nelems > 100 and e%(Nelems//100) == 0:                                     # report in 1% steps if possible and avoid modulo by zero
-    damask.util.progressBar(iteration=e,total=Nelems,start=starttime,prefix='1/3: connecting elements')
+  if options.verbose and Nelems >= 50 and e%(Nelems//50) == 0:                                      # report in 2% steps if possible and avoid modulo by zero
+    damask.util.progressBar(iteration=e,total=Nelems,prefix='1/3: connecting elements')
  for n in map(p.node_sequence,p.element(e).items):
    if n not in elementsOfNode:
      elementsOfNode[n] = [p.element_id(e)]
@ -856,10 +855,9 @@ damask.util.progressBar(iteration=1,total=1,prefix='1/3: connecting elements')

 if options.nodalScalar:
  Npoints = stat['NumberOfNodes']
-  starttime = time.time()
  for n in range(Npoints):
-    if options.verbose and Npoints > 100 and e%(Npoints//100) == 0:                                 # report in 1% steps if possible and avoid modulo by zero
-      damask.util.progressBar(iteration=n,total=Npoints,start=starttime,prefix='2/3: scanning nodes     ')
+    if options.verbose and Npoints >= 50 and e%(Npoints//50) == 0:                                  # report in 2% steps if possible and avoid modulo by zero
+      damask.util.progressBar(iteration=n,total=Npoints,prefix='2/3: scanning nodes     ')
    myNodeID = p.node_id(n)
    myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z]
    myElemID = 0
@ -894,10 +892,9 @@ if options.nodalScalar:

 else:
  Nelems = stat['NumberOfElements']
-  starttime = time.time()
  for e in range(Nelems):
-    if options.verbose and Nelems > 100 and e%(Nelems//100) == 0:                                   # report in 1% steps if possible and avoid modulo by zero
-      damask.util.progressBar(iteration=e,total=Nelems,start=starttime,prefix='2/3: scanning elements  ')
+    if options.verbose and Nelems >= 50 and e%(Nelems//50) == 0:                                    # report in 2% steps if possible and avoid modulo by zero
+      damask.util.progressBar(iteration=e,total=Nelems,prefix='2/3: scanning elements  ')
    myElemID = p.element_id(e)
    myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z],
                                                         list(map(p.node, map(p.node_sequence, p.element(e).items))))))
@ -1035,11 +1032,10 @@ for incCount,position in enumerate(locations):     # walk through locations

  member = 0
  Ngroups = len(groups)
-  starttime = time.time()
  for j,group in enumerate(groups):
    f = incCount*Ngroups + j
-    if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0:                 # report in 1% steps if possible and avoid modulo by zero
-      damask.util.progressBar(iteration=f,total=Ngroups*Nincs,start=starttime,prefix='3/3: processing points  ')
+    if options.verbose and (Ngroups*Nincs) >= 50 and f%((Ngroups*Nincs)//50) == 0:                  # report in 2% steps if possible and avoid modulo by zero
+      damask.util.progressBar(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points  ')
    N = 0                                                                                           # group member counter
    for (e,n,i,g,n_local) in group[1:]:                                                             # loop over group members
      member += 1
@ -1091,7 +1087,7 @@ for incCount,position in enumerate(locations):     # walk through locations
                                      ['Crystallite']*len(options.crystalliteResult) +
                                      ['Constitutive']*len(options.constitutiveResult)
                                      ):
-          outputIndex = (list(zip(*outputFormat[resultType]['outputs']))[0]).index(label)             # find the position of this output in the outputFormat
+          outputIndex = (list(zip(*outputFormat[resultType]['outputs']))[0]).index(label)           # find the position of this output in the outputFormat
          length = int(outputFormat[resultType]['outputs'][outputIndex][1])
          thisHead = heading('_',[[component,''.join( label.split() )] for component in range(int(length>1),length+int(length>1))])
          if assembleHeader: header += thisHead
--- a/python/damask/.gitignore
+++ b/python/damask/.gitignore
--- a/python/damask/init.py
+++ b/python/damask/init.py
--- a/python/damask/asciitable.py
+++ b/python/damask/asciitable.py
--- a/python/damask/colormaps.py
+++ b/python/damask/colormaps.py
--- a/python/damask/config/init.py
+++ b/python/damask/config/init.py
--- a/python/damask/config/material.py
+++ b/python/damask/config/material.py
--- a/python/damask/environment.py
+++ b/python/damask/environment.py
--- a/python/damask/geometry/init.py
+++ b/python/damask/geometry/init.py
--- a/python/damask/geometry/geometry.py
+++ b/python/damask/geometry/geometry.py
--- a/python/damask/geometry/marc.py
+++ b/python/damask/geometry/marc.py
--- a/python/damask/geometry/spectral.py
+++ b/python/damask/geometry/spectral.py
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
--- a/python/damask/result.py
+++ b/python/damask/result.py
--- a/python/damask/result/marc2vtk.py
+++ b/python/damask/result/marc2vtk.py
--- a/python/damask/solver/init.py
+++ b/python/damask/solver/init.py
--- a/python/damask/solver/abaqus.py
+++ b/python/damask/solver/abaqus.py
--- a/python/damask/solver/marc.py
+++ b/python/damask/solver/marc.py
--- a/python/damask/solver/solver.py
+++ b/python/damask/solver/solver.py
--- a/python/damask/solver/spectral.py
+++ b/python/damask/solver/spectral.py
--- a/python/damask/test/init.py
+++ b/python/damask/test/init.py
--- a/python/damask/test/test.py
+++ b/python/damask/test/test.py
--- a/python/damask/util.py
+++ b/python/damask/util.py
@ -134,7 +134,7 @@ class extendableOption(Option):

 # Print iterations progress
 # from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
-def progressBar(iteration, total, start=None, prefix='', suffix='', decimals=1, bar_length=50):
+def progressBar(iteration, total, prefix='', bar_length=50):
  """
  Call in a loop to create terminal progress bar

@ -142,23 +142,29 @@ def progressBar(iteration, total, start=None, prefix='', suffix='', decimals=1,
      iteration   - Required  : current iteration (Int)
      total       - Required  : total iterations (Int)
      prefix      - Optional  : prefix string (Str)
-      suffix      - Optional  : suffix string (Str)
-      decimals    - Optional  : positive number of decimals in percent complete (Int)
      bar_length  - Optional  : character length of bar (Int)
  """
-  import time
-  
-  if suffix == '' and start is not None and iteration > 0:
-    remainder = (total - iteration) * (time.time()-start)/iteration
-    suffix = '{: 3d}:'.format(int( remainder//3600)) + \
-             '{:02d}:'.format(int((remainder//60)%60)) + \
-             '{:02d}' .format(int( remainder     %60))
-  str_format = "{{0:{}.{}f}}".format(decimals+4,decimals)
-  percents = str_format.format(100 * (iteration / float(total)))
-  filled_length = int(round(bar_length * iteration / float(total)))
-  bar = '█' * filled_length + '-' * (bar_length - filled_length)
+  fraction = iteration / float(total)
+  if not hasattr(progressBar, "last_fraction"):                                                     # first call to function
+    progressBar.start_time    = time.time()
+    progressBar.last_fraction = -1.0
+    remaining_time = '   n/a'
+  else:
+    if fraction <= progressBar.last_fraction or iteration == 0:                                     # reset: called within a new loop
+      progressBar.start_time    = time.time()
+      progressBar.last_fraction = -1.0
+      remaining_time = '   n/a'
+    else:
+      progressBar.last_fraction = fraction
+      remainder = (total - iteration) * (time.time()-progressBar.start_time)/iteration
+      remaining_time = '{: 3d}:'.format(int( remainder//3600)) + \
+                       '{:02d}:'.format(int((remainder//60)%60)) + \
+                       '{:02d}' .format(int( remainder     %60))

-  sys.stderr.write('\r%s |%s| %s %s' % (prefix, bar, percents+'%', suffix)),
+  filled_length = int(round(bar_length * fraction))
+  bar = '█' * filled_length + '░' * (bar_length - filled_length)
+
+  sys.stderr.write('\r{} {} {}'.format(prefix, bar, remaining_time)),

  if iteration == total: sys.stderr.write('\n\n')
  sys.stderr.flush()
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@ -90,9 +90,7 @@ list(APPEND OBJECTFILES $<TARGET_OBJECTS:PLASTIC>)
 add_library (KINEMATICS OBJECT
  "kinematics_cleavage_opening.f90"
  "kinematics_slipplane_opening.f90"
-  "kinematics_thermal_expansion.f90"
-  "kinematics_vacancy_strain.f90"
-  "kinematics_hydrogen_strain.f90")
+  "kinematics_thermal_expansion.f90")
 add_dependencies(KINEMATICS DAMASK_HELPERS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:KINEMATICS>)

@ -102,10 +100,7 @@ add_library (SOURCE OBJECT
  "source_damage_isoBrittle.f90"
  "source_damage_isoDuctile.f90"
  "source_damage_anisoBrittle.f90"
-  "source_damage_anisoDuctile.f90"
-  "source_vacancy_phenoplasticity.f90"
-  "source_vacancy_irradiation.f90"
-  "source_vacancy_thermalfluc.f90")
+  "source_damage_anisoDuctile.f90")
 add_dependencies(SOURCE DAMASK_HELPERS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:SOURCE>)

@ -124,25 +119,6 @@ add_library(HOMOGENIZATION OBJECT
 add_dependencies(HOMOGENIZATION CRYSTALLITE)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:HOMOGENIZATION>)

-add_library(HYDROGENFLUX OBJECT
-  "hydrogenflux_isoconc.f90"
-  "hydrogenflux_cahnhilliard.f90")
-add_dependencies(HYDROGENFLUX CRYSTALLITE)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:HYDROGENFLUX>)
-
-add_library(POROSITY OBJECT
-  "porosity_none.f90"
-  "porosity_phasefield.f90")
-add_dependencies(POROSITY CRYSTALLITE)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:POROSITY>)
-
-add_library(VACANCYFLUX OBJECT
-  "vacancyflux_isoconc.f90"
-  "vacancyflux_isochempot.f90"
-  "vacancyflux_cahnhilliard.f90")
-add_dependencies(VACANCYFLUX CRYSTALLITE)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:VACANCYFLUX>)
-
 add_library(DAMAGE OBJECT
  "damage_none.f90"
  "damage_local.f90"
@ -158,7 +134,7 @@ add_dependencies(THERMAL CRYSTALLITE)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:THERMAL>)

 add_library(DAMASK_ENGINE OBJECT "homogenization.f90")
-add_dependencies(DAMASK_ENGINE THERMAL DAMAGE VACANCYFLUX POROSITY HYDROGENFLUX HOMOGENIZATION)
+add_dependencies(DAMASK_ENGINE THERMAL DAMAGE HOMOGENIZATION)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_ENGINE>)

 add_library(DAMASK_CPFE OBJECT "CPFEM2.f90")
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@ -304,8 +304,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
   homogState, &
   thermalState, &
   damageState, &
-   vacancyfluxState, &
-   hydrogenfluxState, &
   phaseAt, phasememberAt, &
   material_phase, &
   phase_plasticity, &
@ -421,8 +419,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
     homogState       (homog)%state0 =  homogState       (homog)%state
     thermalState     (homog)%state0 =  thermalState     (homog)%state
     damageState      (homog)%state0 =  damageState      (homog)%state
-     vacancyfluxState (homog)%state0 =  vacancyfluxState (homog)%state
-     hydrogenfluxState(homog)%state0 =  hydrogenfluxState(homog)%state
   enddo
  

--- a/src/CPFEM2.f90
+++ b/src/CPFEM2.f90
@ -203,8 +203,6 @@ subroutine CPFEM_age()
   homogState, &
   thermalState, &
   damageState, &
-   vacancyfluxState, &
-   hydrogenfluxState, &
   material_phase, &
   phase_plasticity, &
   phase_Nsources
@ -268,8 +266,6 @@ if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
    homogState       (homog)%state0 =  homogState       (homog)%state
    thermalState     (homog)%state0 =  thermalState     (homog)%state
    damageState      (homog)%state0 =  damageState      (homog)%state
-    vacancyfluxState (homog)%state0 =  vacancyfluxState (homog)%state
-    hydrogenfluxState(homog)%state0 =  hydrogenfluxState(homog)%state
  enddo

 if (restartWrite) then
--- a/src/IO.f90
+++ b/src/IO.f90
@ -1251,6 +1251,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
   msg = 'negative number systems requested'
 case (145_pInt)
   msg = 'too many systems requested'
+ case (146_pInt)
+   msg = 'number of values does not match'

 !--------------------------------------------------------------------------------------------------
 ! material error messages and related messages in mesh
--- a/src/commercialFEM_fileList.f90
+++ b/src/commercialFEM_fileList.f90
@ -21,14 +21,9 @@
 #include "source_damage_isoDuctile.f90"
 #include "source_damage_anisoBrittle.f90"
 #include "source_damage_anisoDuctile.f90"
-#include "source_vacancy_phenoplasticity.f90"
-#include "source_vacancy_irradiation.f90"
-#include "source_vacancy_thermalfluc.f90"
 #include "kinematics_cleavage_opening.f90"
 #include "kinematics_slipplane_opening.f90"
 #include "kinematics_thermal_expansion.f90"
-#include "kinematics_vacancy_strain.f90"
-#include "kinematics_hydrogen_strain.f90"
 #include "plastic_none.f90"
 #include "plastic_isotropic.f90"
 #include "plastic_phenopowerlaw.f90"
@ -47,12 +42,5 @@
 #include "damage_none.f90"
 #include "damage_local.f90"
 #include "damage_nonlocal.f90"
-#include "vacancyflux_isoconc.f90"
-#include "vacancyflux_isochempot.f90"
-#include "vacancyflux_cahnhilliard.f90"
-#include "porosity_none.f90"
-#include "porosity_phasefield.f90"
-#include "hydrogenflux_isoconc.f90"
-#include "hydrogenflux_cahnhilliard.f90"
 #include "homogenization.f90"
 #include "CPFEM.f90"
--- a/src/config.f90
+++ b/src/config.f90
@ -550,7 +550,7 @@ end function getString
 !> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
 !! values from the last occurrence. If key is not found exits with error unless default is given.
 !--------------------------------------------------------------------------------------------------
-function getFloats(this,key,defaultVal,requiredShape)
+function getFloats(this,key,defaultVal,requiredShape,requiredSize)
 use IO, only: &
   IO_error, &
   IO_stringValue, &
@ -561,7 +561,8 @@ function getFloats(this,key,defaultVal,requiredShape)
 class(tPartitionedStringList), target, intent(in)   :: this
 character(len=*),              intent(in)           :: key
 real(pReal),   dimension(:),   intent(in), optional :: defaultVal
- integer(pInt), dimension(:),   intent(in), optional :: requiredShape
+ integer(pInt), dimension(:),   intent(in), optional :: requiredShape ! not useful (is always 1D array)
+ integer(pInt),                 intent(in), optional :: requiredSize
 type(tPartitionedStringList),  pointer              :: item
 integer(pInt)                                       :: i
 logical                                             :: found, &
@ -588,6 +589,9 @@ function getFloats(this,key,defaultVal,requiredShape)
 if (.not. found) then
   if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif
 endif
+ if (present(requiredSize)) then
+   if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key)
+ endif

 end function getFloats

@ -597,7 +601,7 @@ end function getFloats
 !> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
 !! values from the last occurrence. If key is not found exits with error unless default is given.
 !--------------------------------------------------------------------------------------------------
-function getInts(this,key,defaultVal,requiredShape)
+function getInts(this,key,defaultVal,requiredShape,requiredSize)
 use IO, only: &
   IO_error, &
   IO_stringValue, &
@ -608,7 +612,8 @@ function getInts(this,key,defaultVal,requiredShape)
 class(tPartitionedStringList), target, intent(in)   :: this
 character(len=*),              intent(in)           :: key
 integer(pInt), dimension(:),   intent(in), optional :: defaultVal, &
-                                                        requiredShape
+                                                        requiredShape ! not useful (is always 1D array)
+ integer(pInt),                 intent(in), optional :: requiredSize
 type(tPartitionedStringList),  pointer              :: item
 integer(pInt)                                       :: i
 logical                                             :: found, &
@ -635,6 +640,9 @@ function getInts(this,key,defaultVal,requiredShape)
 if (.not. found) then
   if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif
 endif
+ if (present(requiredSize)) then
+   if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key)
+ endif

 end function getInts

--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@ -88,14 +88,9 @@ subroutine constitutive_init()
   SOURCE_damage_isoDuctile_ID, &
   SOURCE_damage_anisoBrittle_ID, &
   SOURCE_damage_anisoDuctile_ID, &
-   SOURCE_vacancy_phenoplasticity_ID, &
-   SOURCE_vacancy_irradiation_ID, &
-   SOURCE_vacancy_thermalfluc_ID, &
   KINEMATICS_cleavage_opening_ID, &
   KINEMATICS_slipplane_opening_ID, &
   KINEMATICS_thermal_expansion_ID, &
-   KINEMATICS_vacancy_strain_ID, &
-   KINEMATICS_hydrogen_strain_ID, &
   ELASTICITY_HOOKE_label, &
   PLASTICITY_NONE_label, &
   PLASTICITY_ISOTROPIC_label, &
@ -110,9 +105,6 @@ subroutine constitutive_init()
   SOURCE_damage_isoDuctile_label, &
   SOURCE_damage_anisoBrittle_label, &
   SOURCE_damage_anisoDuctile_label, &
-   SOURCE_vacancy_phenoplasticity_label, &
-   SOURCE_vacancy_irradiation_label, &
-   SOURCE_vacancy_thermalfluc_label, &
   plasticState, &
   sourceState

@ -129,14 +121,9 @@ subroutine constitutive_init()
 use source_damage_isoDuctile
 use source_damage_anisoBrittle
 use source_damage_anisoDuctile
- use source_vacancy_phenoplasticity
- use source_vacancy_irradiation
- use source_vacancy_thermalfluc
 use kinematics_cleavage_opening
 use kinematics_slipplane_opening
 use kinematics_thermal_expansion
- use kinematics_vacancy_strain
- use kinematics_hydrogen_strain

 implicit none
 integer(pInt), parameter :: FILEUNIT = 204_pInt
@ -179,9 +166,6 @@ subroutine constitutive_init()
 if (any(phase_source == SOURCE_damage_isoDuctile_ID))       call source_damage_isoDuctile_init(FILEUNIT)
 if (any(phase_source == SOURCE_damage_anisoBrittle_ID))     call source_damage_anisoBrittle_init(FILEUNIT)
 if (any(phase_source == SOURCE_damage_anisoDuctile_ID))     call source_damage_anisoDuctile_init(FILEUNIT)
- if (any(phase_source == SOURCE_vacancy_phenoplasticity_ID)) call source_vacancy_phenoplasticity_init(FILEUNIT)
- if (any(phase_source == SOURCE_vacancy_irradiation_ID))     call source_vacancy_irradiation_init(FILEUNIT)
- if (any(phase_source == SOURCE_vacancy_thermalfluc_ID))     call source_vacancy_thermalfluc_init(FILEUNIT)

 !--------------------------------------------------------------------------------------------------
 ! parse kinematic mechanisms from config file
@ -189,8 +173,6 @@ subroutine constitutive_init()
 if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID))  call kinematics_cleavage_opening_init(FILEUNIT)
 if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT)
 if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT)
- if (any(phase_kinematics == KINEMATICS_vacancy_strain_ID))    call kinematics_vacancy_strain_init(FILEUNIT)
- if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID))   call kinematics_hydrogen_strain_init(FILEUNIT)
 close(FILEUNIT)

 call config_deallocate('material.config/phase')
@ -283,21 +265,6 @@ subroutine constitutive_init()
             outputName = SOURCE_damage_anisoDuctile_label
             thisOutput => source_damage_anisoDuctile_output
             thisSize   => source_damage_anisoDuctile_sizePostResult
-           case (SOURCE_vacancy_phenoplasticity_ID) sourceType
-             ins = source_vacancy_phenoplasticity_instance(ph)
-             outputName = SOURCE_vacancy_phenoplasticity_label
-             thisOutput => source_vacancy_phenoplasticity_output
-             thisSize   => source_vacancy_phenoplasticity_sizePostResult
-           case (SOURCE_vacancy_irradiation_ID) sourceType
-             ins = source_vacancy_irradiation_instance(ph)
-             outputName = SOURCE_vacancy_irradiation_label
-             thisOutput => source_vacancy_irradiation_output
-             thisSize   => source_vacancy_irradiation_sizePostResult
-           case (SOURCE_vacancy_thermalfluc_ID) sourceType
-             ins = source_vacancy_thermalfluc_instance(ph)
-             outputName = SOURCE_vacancy_thermalfluc_label
-             thisOutput => source_vacancy_thermalfluc_output
-             thisSize   => source_vacancy_thermalfluc_sizePostResult
           case default sourceType
             knownSource = .false.
         end select sourceType
@ -512,8 +479,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
     dLp_dMp = 0.0_pReal

   case (PLASTICITY_ISOTROPIC_ID) plasticityType
-     call plastic_isotropic_LpAndItsTangent       (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
-     dLp_dMp = math_Plain99to3333(dLp_dMp99)                                                        ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_isotropic_LpAndItsTangent       (Lp,dLp_dMp,Mp,instance,of)

   case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
     of = phasememberAt(ipc,ip,el)
@ -560,6 +528,7 @@ end subroutine constitutive_LpAndItsTangents

 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient
+! ToDo: MD: S is Mi?
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el)
 use prec, only: &
@ -568,8 +537,12 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
   math_I3, &
   math_inv33, &
   math_det33, &
-   math_mul33x33
+   math_mul33x33, &
+   math_Mandel6to33
 use material, only: &
+   phasememberAt, &
+   phase_plasticity, &
+   phase_plasticityInstance, &
   phase_plasticity, &
   material_phase, &
   phase_kinematics, &
@ -577,9 +550,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
   PLASTICITY_isotropic_ID, &
   KINEMATICS_cleavage_opening_ID, &
   KINEMATICS_slipplane_opening_ID, &
-   KINEMATICS_thermal_expansion_ID, &
-   KINEMATICS_vacancy_strain_ID, &
-   KINEMATICS_hydrogen_strain_ID
+   KINEMATICS_thermal_expansion_ID
 use plastic_isotropic, only: &
   plastic_isotropic_LiAndItsTangent
 use kinematics_cleavage_opening, only: &
@ -588,10 +559,6 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
   kinematics_slipplane_opening_LiAndItsTangent
 use kinematics_thermal_expansion, only: &
   kinematics_thermal_expansion_LiAndItsTangent
- use kinematics_vacancy_strain, only: &
-   kinematics_vacancy_strain_LiAndItsTangent
- use kinematics_hydrogen_strain, only: &
-   kinematics_hydrogen_strain_LiAndItsTangent

 implicit none
 integer(pInt), intent(in) :: &
@ -607,19 +574,18 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
 real(pReal),   intent(out), dimension(3,3,3,3) :: &
   dLi_dS, &                                                                                        !< derivative of Li with respect to S
   dLi_dFi
+ 
 real(pReal), dimension(3,3) :: &
-   my_Li                                                                                            !< intermediate velocity gradient
- real(pReal), dimension(3,3,3,3) :: &
-   my_dLi_dS
- real(pReal), dimension(3,3) :: &
+   my_Li, &                                                                                            !< intermediate velocity gradient
   FiInv, &
   temp_33
+ real(pReal), dimension(3,3,3,3) :: &
+   my_dLi_dS
 real(pReal) :: &
   detFi
 integer(pInt) :: &
-   k                                                                                                !< counter in kinematics loop
- integer(pInt) :: &
-   i, j
+   k, i, j, &
+   instance, of

 Li = 0.0_pReal
 dLi_dS  = 0.0_pReal
@ -627,7 +593,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e

 plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
   case (PLASTICITY_isotropic_ID) plasticityType
-     call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_Mandel6to33(S6),instance,of)
   case default plasticityType
     my_Li = 0.0_pReal
     my_dLi_dS = 0.0_pReal
@ -644,10 +612,6 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
       call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
     case (KINEMATICS_thermal_expansion_ID) kinematicsType
       call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
-     case (KINEMATICS_vacancy_strain_ID) kinematicsType
-       call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
-     case (KINEMATICS_hydrogen_strain_ID) kinematicsType
-       call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
     case default kinematicsType
       my_Li = 0.0_pReal
       my_dLi_dS = 0.0_pReal
@ -684,15 +648,9 @@ pure function constitutive_initialFi(ipc, ip, el)
   phase_kinematics, &
   phase_Nkinematics, &
   material_phase, &
-   KINEMATICS_thermal_expansion_ID, &
-   KINEMATICS_vacancy_strain_ID, &
-   KINEMATICS_hydrogen_strain_ID
+   KINEMATICS_thermal_expansion_ID
 use kinematics_thermal_expansion, only: &
   kinematics_thermal_expansion_initialStrain
- use kinematics_vacancy_strain, only: &
-   kinematics_vacancy_strain_initialStrain
- use kinematics_hydrogen_strain, only: &
-   kinematics_hydrogen_strain_initialStrain

 implicit none
 integer(pInt), intent(in) :: &
@ -711,12 +669,6 @@ pure function constitutive_initialFi(ipc, ip, el)
     case (KINEMATICS_thermal_expansion_ID) kinematicsType
       constitutive_initialFi = &
         constitutive_initialFi + kinematics_thermal_expansion_initialStrain(ipc, ip, el)
-     case (KINEMATICS_vacancy_strain_ID) kinematicsType
-       constitutive_initialFi = &
-         constitutive_initialFi + kinematics_vacancy_strain_initialStrain(ipc, ip, el)
-     case (KINEMATICS_hydrogen_strain_ID) kinematicsType
-       constitutive_initialFi = &
-         constitutive_initialFi + kinematics_hydrogen_strain_initialStrain(ipc, ip, el)
   end select kinematicsType
 enddo KinematicsLoop

@ -771,10 +723,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
   phase_stiffnessDegradation, &
   damage, &
   damageMapping, &
-   porosity, &
-   porosityMapping, &
-   STIFFNESS_DEGRADATION_damage_ID, &
-   STIFFNESS_DEGRADATION_porosity_ID
+   STIFFNESS_DEGRADATION_damage_ID

 implicit none
 integer(pInt), intent(in) :: &
@ -804,8 +753,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
   degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el)))
     case (STIFFNESS_DEGRADATION_damage_ID) degradationType
       C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2_pInt
-     case (STIFFNESS_DEGRADATION_porosity_ID) degradationType
-       C = C * porosity(ho)%p(porosityMapping(ho)%p(ip,el))**2_pInt
   end select degradationType
 enddo DegradationLoop

@ -916,7 +863,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
 plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))

   case (PLASTICITY_ISOTROPIC_ID) plasticityType
-     call plastic_isotropic_dotState    (math_Mandel33to6(Mp),ipc,ip,el)
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_isotropic_dotState    (Mp,instance,of)

   case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
     of = phasememberAt(ipc,ip,el)
@ -986,19 +935,13 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
   material_phase, &
   PLASTICITY_KINEHARDENING_ID, &
   PLASTICITY_NONLOCAL_ID, &
-   SOURCE_damage_isoBrittle_ID, &
-   SOURCE_vacancy_irradiation_ID, &
-   SOURCE_vacancy_thermalfluc_ID
+   SOURCE_damage_isoBrittle_ID
 use plastic_kinehardening, only: &
   plastic_kinehardening_deltaState   
 use plastic_nonlocal, only: &
   plastic_nonlocal_deltaState
 use source_damage_isoBrittle, only: &
   source_damage_isoBrittle_deltaState
- use source_vacancy_irradiation, only: &
-   source_vacancy_irradiation_deltaState
- use source_vacancy_thermalfluc, only: &
-   source_vacancy_thermalfluc_deltaState

 implicit none
 integer(pInt), intent(in) :: &
@ -1035,12 +978,6 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
       call source_damage_isoBrittle_deltaState  (constitutive_homogenizedC(ipc,ip,el), Fe, &
                                                  ipc, ip, el)

-     case (SOURCE_vacancy_irradiation_ID) sourceType
-       call source_vacancy_irradiation_deltaState(ipc, ip, el)
-
-     case (SOURCE_vacancy_thermalfluc_ID) sourceType
-       call source_vacancy_thermalfluc_deltaState(ipc, ip, el)
-
   end select sourceType

 enddo SourceLoop
@ -1140,8 +1077,10 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)

 plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
   case (PLASTICITY_ISOTROPIC_ID) plasticityType
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
     constitutive_postResults(startPos:endPos) = &
-       plastic_isotropic_postResults(S6,ipc,ip,el)
+       plastic_isotropic_postResults(Mp,instance,of)

   case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
     of = phasememberAt(ipc,ip,el)
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@ -25,10 +25,7 @@ module homogenization
   materialpoint_sizeResults, &
   homogenization_maxSizePostResults, &
   thermal_maxSizePostResults, &
-   damage_maxSizePostResults, &
-   vacancyflux_maxSizePostResults, &
-   porosity_maxSizePostResults, &
-   hydrogenflux_maxSizePostResults
+   damage_maxSizePostResults

 real(pReal),   dimension(:,:,:,:),     allocatable, private :: &
   materialpoint_subF0, &                                                                           !< def grad of IP at beginning of homogenization increment
@ -100,13 +97,6 @@ subroutine homogenization_init
 use damage_none
 use damage_local
 use damage_nonlocal
- use vacancyflux_isoconc
- use vacancyflux_isochempot
- use vacancyflux_cahnhilliard
- use porosity_none
- use porosity_phasefield
- use hydrogenflux_isoconc
- use hydrogenflux_cahnhilliard
 use IO
 use numerics, only: &
   worldrank
@ -155,33 +145,6 @@ subroutine homogenization_init
 if (any(damage_type == DAMAGE_nonlocal_ID)) &
   call damage_nonlocal_init(FILEUNIT)

-!--------------------------------------------------------------------------------------------------
-! parse vacancy transport from config file
- call IO_checkAndRewind(FILEUNIT)
- if (any(vacancyflux_type == VACANCYFLUX_isoconc_ID)) &
-   call vacancyflux_isoconc_init()
- if (any(vacancyflux_type == VACANCYFLUX_isochempot_ID)) &
-   call vacancyflux_isochempot_init(FILEUNIT)
- if (any(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID)) &
-   call vacancyflux_cahnhilliard_init(FILEUNIT)
-
-!--------------------------------------------------------------------------------------------------
-! parse porosity from config file
- call IO_checkAndRewind(FILEUNIT)
- if (any(porosity_type == POROSITY_none_ID)) &
-   call porosity_none_init()
- if (any(porosity_type == POROSITY_phasefield_ID)) &
-   call porosity_phasefield_init(FILEUNIT)
-
-!--------------------------------------------------------------------------------------------------
-! parse hydrogen transport from config file
- call IO_checkAndRewind(FILEUNIT)
- if (any(hydrogenflux_type == HYDROGENFLUX_isoconc_ID)) &
-   call hydrogenflux_isoconc_init()
- if (any(hydrogenflux_type == HYDROGENFLUX_cahnhilliard_ID)) &
-   call hydrogenflux_cahnhilliard_init(FILEUNIT)
- close(FILEUNIT)
-
 !--------------------------------------------------------------------------------------------------
 ! write description file for homogenization output
 mainProcess2: if (worldrank == 0) then
@ -277,83 +240,6 @@ subroutine homogenization_init
           enddo
         endif
       endif
-       i = vacancyflux_typeInstance(p)                                                                  ! which instance of this vacancy flux type
-       valid = .true.                                                                                   ! assume valid
-       select case(vacancyflux_type(p))                                                                 ! split per vacancy flux type
-         case (VACANCYFLUX_isoconc_ID)
-           outputName = VACANCYFLUX_isoconc_label
-           thisNoutput => null()
-           thisOutput => null()
-           thisSize   => null()
-         case (VACANCYFLUX_isochempot_ID)
-           outputName = VACANCYFLUX_isochempot_label
-           thisNoutput => vacancyflux_isochempot_Noutput
-           thisOutput => vacancyflux_isochempot_output
-           thisSize   => vacancyflux_isochempot_sizePostResult
-         case (VACANCYFLUX_cahnhilliard_ID)
-           outputName = VACANCYFLUX_cahnhilliard_label
-           thisNoutput => vacancyflux_cahnhilliard_Noutput
-           thisOutput => vacancyflux_cahnhilliard_output
-           thisSize   => vacancyflux_cahnhilliard_sizePostResult
-         case default
-           valid = .false.
-       end select
-       if (valid) then
-         write(FILEUNIT,'(a)') '(vacancyflux)'//char(9)//trim(outputName)
-         if (vacancyflux_type(p) /= VACANCYFLUX_isoconc_ID) then
-           do e = 1,thisNoutput(i)
-             write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
-           enddo
-         endif
-       endif
-       i = porosity_typeInstance(p)                                                                     ! which instance of this porosity type
-       valid = .true.                                                                                   ! assume valid
-       select case(porosity_type(p))                                                                    ! split per porosity type
-         case (POROSITY_none_ID)
-           outputName = POROSITY_none_label
-           thisNoutput => null()
-           thisOutput => null()
-           thisSize   => null()
-         case (POROSITY_phasefield_ID)
-           outputName = POROSITY_phasefield_label
-           thisNoutput => porosity_phasefield_Noutput
-           thisOutput => porosity_phasefield_output
-           thisSize   => porosity_phasefield_sizePostResult
-         case default
-           valid = .false.
-       end select
-       if (valid) then
-         write(FILEUNIT,'(a)') '(porosity)'//char(9)//trim(outputName)
-         if (porosity_type(p) /= POROSITY_none_ID) then
-           do e = 1,thisNoutput(i)
-             write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
-           enddo
-         endif
-       endif
-       i = hydrogenflux_typeInstance(p)                                                                 ! which instance of this hydrogen flux type
-       valid = .true.                                                                                   ! assume valid
-       select case(hydrogenflux_type(p))                                                                ! split per hydrogen flux type
-         case (HYDROGENFLUX_isoconc_ID)
-           outputName = HYDROGENFLUX_isoconc_label
-           thisNoutput => null()
-           thisOutput => null()
-           thisSize   => null()
-         case (HYDROGENFLUX_cahnhilliard_ID)
-           outputName = HYDROGENFLUX_cahnhilliard_label
-           thisNoutput => hydrogenflux_cahnhilliard_Noutput
-           thisOutput => hydrogenflux_cahnhilliard_output
-           thisSize   => hydrogenflux_cahnhilliard_sizePostResult
-         case default
-           valid = .false.
-       end select
-       if (valid) then
-         write(FILEUNIT,'(a)') '(hydrogenflux)'//char(9)//trim(outputName)
-         if (hydrogenflux_type(p) /= HYDROGENFLUX_isoconc_ID) then
-           do e = 1,thisNoutput(i)
-             write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
-           enddo
-         endif
-       endif
     endif
   enddo
   close(FILEUNIT)
@ -383,25 +269,16 @@ subroutine homogenization_init
 homogenization_maxSizePostResults = 0_pInt
 thermal_maxSizePostResults        = 0_pInt
 damage_maxSizePostResults         = 0_pInt
- vacancyflux_maxSizePostResults    = 0_pInt
- porosity_maxSizePostResults       = 0_pInt
- hydrogenflux_maxSizePostResults   = 0_pInt
 do p = 1,size(config_homogenization)
   homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState       (p)%sizePostResults)
   thermal_maxSizePostResults        = max(thermal_maxSizePostResults,       thermalState     (p)%sizePostResults)
   damage_maxSizePostResults         = max(damage_maxSizePostResults        ,damageState      (p)%sizePostResults)
-   vacancyflux_maxSizePostResults    = max(vacancyflux_maxSizePostResults   ,vacancyfluxState (p)%sizePostResults)
-   porosity_maxSizePostResults       = max(porosity_maxSizePostResults      ,porosityState    (p)%sizePostResults)
-   hydrogenflux_maxSizePostResults   = max(hydrogenflux_maxSizePostResults  ,hydrogenfluxState(p)%sizePostResults)
 enddo

 materialpoint_sizeResults = 1 &                                                                    ! grain count
                           + 1 + homogenization_maxSizePostResults &                                ! homogSize & homogResult
                               + thermal_maxSizePostResults        &
                               + damage_maxSizePostResults         &
-                               + vacancyflux_maxSizePostResults    &
-                               + porosity_maxSizePostResults       &
-                               + hydrogenflux_maxSizePostResults   &
                           + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults &      ! crystallite size & crystallite results
                                                        + 1 + constitutive_plasticity_maxSizePostResults &     ! constitutive size & constitutive results
                                                            + constitutive_source_maxSizePostResults)
@ -460,9 +337,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
   homogState, &
   thermalState, &
   damageState, &
-   vacancyfluxState, &
-   porosityState, &
-   hydrogenfluxState, &
   phase_Nsources, &
   mappingHomogenization, &
   phaseAt, phasememberAt, &
@ -569,18 +443,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
     damageState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
       damageState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
       damageState(mappingHomogenization(2,i,e))%State0(   :,mappingHomogenization(1,i,e))          ! ...internal damage state
-   forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-     vacancyfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-       vacancyfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
-       vacancyfluxState(mappingHomogenization(2,i,e))%State0(   :,mappingHomogenization(1,i,e))     ! ...internal vacancy transport state
-   forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-     porosityState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-       porosityState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
-       porosityState(mappingHomogenization(2,i,e))%State0(   :,mappingHomogenization(1,i,e))        ! ...internal porosity state
-   forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-     hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-       hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
-       hydrogenfluxState(mappingHomogenization(2,i,e))%State0(   :,mappingHomogenization(1,i,e))    ! ...internal hydrogen transport state
 enddo
 NiterationHomog = 0_pInt

@ -654,18 +516,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
             damageState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
               damageState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
               damageState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e))  ! ...internal damage state
-           forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-             vacancyfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-               vacancyfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
-               vacancyfluxState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e))! ...internal vacancy transport state
-           forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-             porosityState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-               porosityState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
-               porosityState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e))! ...internal porosity state
-           forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-             hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-               hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
-               hydrogenfluxState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state
           materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e)                       ! ...def grad
         endif steppingNeeded

@ -729,18 +579,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
             damageState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
               damageState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e)) = &
               damageState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))  ! ...internal damage state
-           forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-             vacancyfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-               vacancyfluxState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e)) = &
-               vacancyfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal vacancy transport state
-           forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-             porosityState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-               porosityState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e)) = &
-               porosityState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal porosity state
-           forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
-             hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
-               hydrogenfluxState(mappingHomogenization(2,i,e))%State(    :,mappingHomogenization(1,i,e)) = &
-               hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal hydrogen transport state
         endif
       endif converged

@ -846,9 +684,6 @@ subroutine materialpoint_postResults
   homogState, &
   thermalState, &
   damageState, &
-   vacancyfluxState, &
-   porosityState, &
-   hydrogenfluxState, &
   plasticState, &
   sourceState, &
   material_phase, &
@ -877,10 +712,7 @@ subroutine materialpoint_postResults

       theSize = homogState       (mappingHomogenization(2,i,e))%sizePostResults &
               + thermalState     (mappingHomogenization(2,i,e))%sizePostResults &
-               + damageState      (mappingHomogenization(2,i,e))%sizePostResults &
-               + vacancyfluxState (mappingHomogenization(2,i,e))%sizePostResults &
-               + porosityState    (mappingHomogenization(2,i,e))%sizePostResults &
-               + hydrogenfluxState(mappingHomogenization(2,i,e))%sizePostResults
+               + damageState      (mappingHomogenization(2,i,e))%sizePostResults
       materialpoint_results(thePos+1,i,e) = real(theSize,pReal)                                    ! tell size of homogenization results
       thePos = thePos + 1_pInt

@ -964,12 +796,10 @@ function homogenization_updateState(ip,el)
   homogenization_type, &
   thermal_type, &
   damage_type, &
-   vacancyflux_type, &
   homogenization_maxNgrains, &
   HOMOGENIZATION_RGC_ID, &
   THERMAL_adiabatic_ID, &
-   DAMAGE_local_ID, &
-   VACANCYFLUX_isochempot_ID
+   DAMAGE_local_ID
 use crystallite, only: &
   crystallite_P, &
   crystallite_dPdF, &
@ -981,8 +811,6 @@ function homogenization_updateState(ip,el)
   thermal_adiabatic_updateState
 use damage_local, only: &
   damage_local_updateState
- use vacancyflux_isochempot, only: &
-   vacancyflux_isochempot_updateState

 implicit none
 integer(pInt), intent(in) :: &
@ -1023,15 +851,6 @@ function homogenization_updateState(ip,el)
                                el)
 end select chosenDamage

- chosenVacancyflux: select case (vacancyflux_type(mesh_element(3,el)))
-   case (VACANCYFLUX_isochempot_ID) chosenVacancyflux
-     homogenization_updateState = &
-       homogenization_updateState .and. &
-       vacancyflux_isochempot_updateState(materialpoint_subdt(ip,el), &
-                                          ip, &
-                                          el)
- end select chosenVacancyflux
-
 end function homogenization_updateState


@ -1095,15 +914,9 @@ function homogenization_postResults(ip,el)
   homogState, &
   thermalState, &
   damageState, &
-   vacancyfluxState, &
-   porosityState, &
-   hydrogenfluxState, &
   homogenization_type, &
   thermal_type, &
   damage_type, &
-   vacancyflux_type, &
-   porosity_type, &
-   hydrogenflux_type, &
   HOMOGENIZATION_NONE_ID, &
   HOMOGENIZATION_ISOSTRAIN_ID, &
   HOMOGENIZATION_RGC_ID, &
@ -1112,14 +925,7 @@ function homogenization_postResults(ip,el)
   THERMAL_conduction_ID, &
   DAMAGE_none_ID, &
   DAMAGE_local_ID, &
-   DAMAGE_nonlocal_ID, &
-   VACANCYFLUX_isoconc_ID, &
-   VACANCYFLUX_isochempot_ID, &
-   VACANCYFLUX_cahnhilliard_ID, &
-   POROSITY_none_ID, &
-   POROSITY_phasefield_ID, &
-   HYDROGENFLUX_isoconc_ID, &
-   HYDROGENFLUX_cahnhilliard_ID
+   DAMAGE_nonlocal_ID
 use homogenization_isostrain, only: &
   homogenization_isostrain_postResults
 use homogenization_RGC, only: &
@ -1132,14 +938,6 @@ function homogenization_postResults(ip,el)
   damage_local_postResults
 use damage_nonlocal, only: &
   damage_nonlocal_postResults
- use vacancyflux_isochempot, only: &
-   vacancyflux_isochempot_postResults
- use vacancyflux_cahnhilliard, only: &
-   vacancyflux_cahnhilliard_postResults
- use porosity_phasefield, only: &
-   porosity_phasefield_postResults
- use hydrogenflux_cahnhilliard, only: &
-   hydrogenflux_cahnhilliard_postResults

 implicit none
 integer(pInt), intent(in) :: &
@ -1147,10 +945,7 @@ function homogenization_postResults(ip,el)
   el                                                                                               !< element number
 real(pReal), dimension(  homogState       (mappingHomogenization(2,ip,el))%sizePostResults &
                        + thermalState     (mappingHomogenization(2,ip,el))%sizePostResults &
-                        + damageState      (mappingHomogenization(2,ip,el))%sizePostResults &
-                        + vacancyfluxState (mappingHomogenization(2,ip,el))%sizePostResults &
-                        + porosityState    (mappingHomogenization(2,ip,el))%sizePostResults &
-                        + hydrogenfluxState(mappingHomogenization(2,ip,el))%sizePostResults) :: &
+                        + damageState      (mappingHomogenization(2,ip,el))%sizePostResults) :: &
   homogenization_postResults
 integer(pInt) :: &
   startPos, endPos
@ -1205,39 +1000,6 @@ function homogenization_postResults(ip,el)
       damage_nonlocal_postResults(ip, el)
 end select chosenDamage

- startPos = endPos + 1_pInt
- endPos   = endPos + vacancyfluxState(mappingHomogenization(2,ip,el))%sizePostResults
- chosenVacancyflux: select case (vacancyflux_type(mesh_element(3,el)))
-   case (VACANCYFLUX_isoconc_ID) chosenVacancyflux
-
-   case (VACANCYFLUX_isochempot_ID) chosenVacancyflux
-     homogenization_postResults(startPos:endPos) = &
-       vacancyflux_isochempot_postResults(ip, el)
-   case (VACANCYFLUX_cahnhilliard_ID) chosenVacancyflux
-     homogenization_postResults(startPos:endPos) = &
-       vacancyflux_cahnhilliard_postResults(ip, el)
- end select chosenVacancyflux
-
- startPos = endPos + 1_pInt
- endPos   = endPos + porosityState(mappingHomogenization(2,ip,el))%sizePostResults
- chosenPorosity: select case (porosity_type(mesh_element(3,el)))
-   case (POROSITY_none_ID) chosenPorosity
-
-   case (POROSITY_phasefield_ID) chosenPorosity
-     homogenization_postResults(startPos:endPos) = &
-       porosity_phasefield_postResults(ip, el)
- end select chosenPorosity
-
- startPos = endPos + 1_pInt
- endPos   = endPos + hydrogenfluxState(mappingHomogenization(2,ip,el))%sizePostResults
- chosenHydrogenflux: select case (hydrogenflux_type(mesh_element(3,el)))
-   case (HYDROGENFLUX_isoconc_ID) chosenHydrogenflux
-
-   case (HYDROGENFLUX_cahnhilliard_ID) chosenHydrogenflux
-     homogenization_postResults(startPos:endPos) = &
-       hydrogenflux_cahnhilliard_postResults(ip, el)
- end select chosenHydrogenflux
-
 end function homogenization_postResults

 end module homogenization
--- a/src/hydrogenflux_cahnhilliard.f90
+++ b/src/hydrogenflux_cahnhilliard.f90
@ -1,508 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for conservative transport of solute hydrogen
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module hydrogenflux_cahnhilliard
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   hydrogenflux_cahnhilliard_sizePostResults                                                           !< cumulative size of post results
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   hydrogenflux_cahnhilliard_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   hydrogenflux_cahnhilliard_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   hydrogenflux_cahnhilliard_Noutput                                                                   !< number of outputs per instance of this damage 
-
- real(pReal),                                                 parameter,           private :: &
-   kB = 1.3806488e-23_pReal                                                                          !< Boltzmann constant in J/Kelvin
-
- enum, bind(c) 
-   enumerator :: undefined_ID, &
-                 hydrogenConc_ID
- end enum
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
-   hydrogenflux_cahnhilliard_outputID                                                                  !< ID of each post result output
-
-
- public :: &
-   hydrogenflux_cahnhilliard_init, &
-   hydrogenflux_cahnhilliard_getMobility33, &
-   hydrogenflux_cahnhilliard_getDiffusion33, &
-   hydrogenflux_cahnhilliard_getFormationEnergy, &
-   hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent, &
-   hydrogenflux_cahnhilliard_getChemPotAndItsTangent, &
-   hydrogenflux_cahnhilliard_putHydrogenConcAndItsRate, &
-   hydrogenflux_cahnhilliard_postResults
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine hydrogenflux_cahnhilliard_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   hydrogenflux_type, &
-   hydrogenflux_typeInstance, &
-   homogenization_Noutput, &
-   HYDROGENFLUX_cahnhilliard_label, &
-   HYDROGENFLUX_cahnhilliard_ID, &
-   material_homog, &  
-   mappingHomogenization, &
-   hydrogenfluxState, &
-   hydrogenfluxMapping, &
-   hydrogenConc, &
-   hydrogenConcRate, &
-   hydrogenflux_initialCh
- use config, only: &
-   material_partHomogenization, &
-   material_partPhase
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o
- integer(pInt) :: sizeState
- integer(pInt) :: NofMyHomog   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  hydrogenflux_'//HYDROGENFLUX_cahnhilliard_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(hydrogenflux_type == HYDROGENFLUX_cahnhilliard_ID),pInt)
- if (maxNinstance == 0_pInt) return
- 
- allocate(hydrogenflux_cahnhilliard_sizePostResults(maxNinstance),                               source=0_pInt)
- allocate(hydrogenflux_cahnhilliard_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
- allocate(hydrogenflux_cahnhilliard_output         (maxval(homogenization_Noutput),maxNinstance))
-          hydrogenflux_cahnhilliard_output = ''
- allocate(hydrogenflux_cahnhilliard_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
- allocate(hydrogenflux_cahnhilliard_Noutput        (maxNinstance),                               source=0_pInt) 
-
- rewind(fileUnit)
- section = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingHomog: do while (trim(line) /= IO_EOF)                                                      ! read through sections of homog part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next homog section
-     section = section + 1_pInt                                                                     ! advance homog section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (section > 0_pInt ) then; if (hydrogenflux_type(section) == HYDROGENFLUX_cahnhilliard_ID) then  ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = hydrogenflux_typeInstance(section)                                                   ! which instance of my hydrogenflux is present homog
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('hydrogenconc')
-             hydrogenflux_cahnhilliard_Noutput(instance) = hydrogenflux_cahnhilliard_Noutput(instance) + 1_pInt
-             hydrogenflux_cahnhilliard_outputID(hydrogenflux_cahnhilliard_Noutput(instance),instance) = hydrogenConc_ID
-             hydrogenflux_cahnhilliard_output(hydrogenflux_cahnhilliard_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
-
-     end select
-   endif; endif
- enddo parsingHomog
- 
- rewind(fileUnit)
- section = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase)         ! wind forward to <homogenization>
-   line = IO_read(fileUnit)
- enddo
- 
- initializeInstances: do section = 1_pInt, size(hydrogenflux_type)
-   if (hydrogenflux_type(section) == HYDROGENFLUX_cahnhilliard_ID) then
-     NofMyHomog=count(material_homog==section)
-     instance = hydrogenflux_typeInstance(section)
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,hydrogenflux_cahnhilliard_Noutput(instance)
-       select case(hydrogenflux_cahnhilliard_outputID(o,instance))
-         case(hydrogenConc_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          hydrogenflux_cahnhilliard_sizePostResult(o,instance) = mySize
-          hydrogenflux_cahnhilliard_sizePostResults(instance)  = hydrogenflux_cahnhilliard_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-
-! allocate state arrays
-     sizeState = 0_pInt
-     hydrogenfluxState(section)%sizeState = sizeState
-     hydrogenfluxState(section)%sizePostResults = hydrogenflux_cahnhilliard_sizePostResults(instance)
-     allocate(hydrogenfluxState(section)%state0   (sizeState,NofMyHomog))
-     allocate(hydrogenfluxState(section)%subState0(sizeState,NofMyHomog))
-     allocate(hydrogenfluxState(section)%state    (sizeState,NofMyHomog))
-
-     nullify(hydrogenfluxMapping(section)%p)
-     hydrogenfluxMapping(section)%p => mappingHomogenization(1,:,:)
-     deallocate(hydrogenConc    (section)%p)
-     deallocate(hydrogenConcRate(section)%p)
-     allocate  (hydrogenConc    (section)%p(NofMyHomog), source=hydrogenflux_initialCh(section))
-     allocate  (hydrogenConcRate(section)%p(NofMyHomog), source=0.0_pReal)
-     
-   endif
- 
- enddo initializeInstances
- 
-end subroutine hydrogenflux_cahnhilliard_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized solute mobility tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function hydrogenflux_cahnhilliard_getMobility33(ip,el)
- use lattice, only: &
-   lattice_hydrogenfluxMobility33
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_push33ToRef
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   hydrogenflux_cahnhilliard_getMobility33
- integer(pInt) :: &
-   grain
-  
- hydrogenflux_cahnhilliard_getMobility33 = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   hydrogenflux_cahnhilliard_getMobility33 = hydrogenflux_cahnhilliard_getMobility33 + &
-    crystallite_push33ToRef(grain,ip,el,lattice_hydrogenfluxMobility33(:,:,material_phase(grain,ip,el)))
- enddo
-
- hydrogenflux_cahnhilliard_getMobility33 = &
-   hydrogenflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function hydrogenflux_cahnhilliard_getMobility33
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized solute nonlocal diffusion tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function hydrogenflux_cahnhilliard_getDiffusion33(ip,el)
- use lattice, only: &
-   lattice_hydrogenfluxDiffusion33
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_push33ToRef
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   hydrogenflux_cahnhilliard_getDiffusion33
- integer(pInt) :: &
-   grain
-  
- hydrogenflux_cahnhilliard_getDiffusion33 = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   hydrogenflux_cahnhilliard_getDiffusion33 = hydrogenflux_cahnhilliard_getDiffusion33 + &
-    crystallite_push33ToRef(grain,ip,el,lattice_hydrogenfluxDiffusion33(:,:,material_phase(grain,ip,el)))
- enddo
-
- hydrogenflux_cahnhilliard_getDiffusion33 = &
-   hydrogenflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function hydrogenflux_cahnhilliard_getDiffusion33
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized solution energy
-!--------------------------------------------------------------------------------------------------
-function hydrogenflux_cahnhilliard_getFormationEnergy(ip,el)
- use lattice, only: &
-   lattice_hydrogenFormationEnergy, &
-   lattice_hydrogenVol, &
-   lattice_hydrogenSurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal) :: &
-   hydrogenflux_cahnhilliard_getFormationEnergy
- integer(pInt) :: &
-   grain
-  
- hydrogenflux_cahnhilliard_getFormationEnergy = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   hydrogenflux_cahnhilliard_getFormationEnergy = hydrogenflux_cahnhilliard_getFormationEnergy + &
-    lattice_hydrogenFormationEnergy(material_phase(grain,ip,el))/ &
-    lattice_hydrogenVol(material_phase(grain,ip,el))/ &
-    lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el))
- enddo
-
- hydrogenflux_cahnhilliard_getFormationEnergy = &
-   hydrogenflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function hydrogenflux_cahnhilliard_getFormationEnergy
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized hydrogen entropy coefficient
-!--------------------------------------------------------------------------------------------------
-function hydrogenflux_cahnhilliard_getEntropicCoeff(ip,el)
- use lattice, only: &
-   lattice_hydrogenVol, &
-   lattice_hydrogenSurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_homog, &
-   material_phase, &
-   temperature, &
-   thermalMapping
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal) :: &
-   hydrogenflux_cahnhilliard_getEntropicCoeff
- integer(pInt) :: &
-   grain
-  
- hydrogenflux_cahnhilliard_getEntropicCoeff = 0.0_pReal
- do grain = 1, homogenization_Ngrains(material_homog(ip,el))
-   hydrogenflux_cahnhilliard_getEntropicCoeff = hydrogenflux_cahnhilliard_getEntropicCoeff + &
-    kB/ &
-    lattice_hydrogenVol(material_phase(grain,ip,el))/ &
-    lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el))
- enddo
-
- hydrogenflux_cahnhilliard_getEntropicCoeff = hydrogenflux_cahnhilliard_getEntropicCoeff* &
-   temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ &
-   real(homogenization_Ngrains(material_homog(ip,el)),pReal)
- 
-end function hydrogenflux_cahnhilliard_getEntropicCoeff
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized kinematic contribution to chemical potential
-!--------------------------------------------------------------------------------------------------
-subroutine hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCh, Ch, ip, el)
- use lattice, only: &
-   lattice_hydrogenSurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_homog, &
-   phase_kinematics, &
-   phase_Nkinematics, &
-   material_phase, &
-   KINEMATICS_hydrogen_strain_ID
- use crystallite, only: &
-   crystallite_Tstar_v, &
-   crystallite_Fi0, &
-   crystallite_Fi
- use kinematics_hydrogen_strain, only: &
-   kinematics_hydrogen_strain_ChemPotAndItsTangent
-
- implicit none
- integer(pInt), intent(in)  :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in)  :: &
-   Ch
- real(pReal),   intent(out) :: &
-   KPot, dKPot_dCh
- real(pReal) :: &
-   my_KPot, my_dKPot_dCh
- integer(pInt) :: &
-   grain, kinematics
-  
- KPot = 0.0_pReal
- dKPot_dCh = 0.0_pReal
- do grain = 1_pInt,homogenization_Ngrains(material_homog(ip,el)) 
-   do kinematics = 1_pInt, phase_Nkinematics(material_phase(grain,ip,el))
-     select case (phase_kinematics(kinematics,material_phase(grain,ip,el)))
-       case (KINEMATICS_hydrogen_strain_ID)
-         call kinematics_hydrogen_strain_ChemPotAndItsTangent(my_KPot, my_dKPot_dCh, &
-                                                             crystallite_Tstar_v(1:6,grain,ip,el), &
-                                                             crystallite_Fi0(1:3,1:3,grain,ip,el), &
-                                                             crystallite_Fi (1:3,1:3,grain,ip,el), &
-                                                             grain,ip, el)
-
-       case default
-         my_KPot = 0.0_pReal
-         my_dKPot_dCh = 0.0_pReal
-     
-     end select
-     KPot = KPot + my_KPot/lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el))
-     dKPot_dCh = dKPot_dCh + my_dKPot_dCh/lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el))
-   enddo
- enddo 
- 
- KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
- dKPot_dCh = dKPot_dCh/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
-
-end subroutine hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized chemical potential
-!--------------------------------------------------------------------------------------------------
-subroutine hydrogenflux_cahnhilliard_getChemPotAndItsTangent(ChemPot,dChemPot_dCh,Ch,ip,el)
- use numerics, only: &
-   hydrogenBoundPenalty, &
-   hydrogenPolyOrder
-
- implicit none
- integer(pInt), intent(in)  :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in)  :: &
-   Ch
- real(pReal),   intent(out) :: &
-   ChemPot, &
-   dChemPot_dCh
- real(pReal) :: &
-   kBT, KPot, dKPot_dCh  
- integer(pInt) :: &
-   o  
-  
- ChemPot = hydrogenflux_cahnhilliard_getFormationEnergy(ip,el)
- dChemPot_dCh = 0.0_pReal
- kBT = hydrogenflux_cahnhilliard_getEntropicCoeff(ip,el)
- do o = 1_pInt, hydrogenPolyOrder
-   ChemPot = ChemPot + kBT*((2.0_pReal*Ch - 1.0_pReal)**real(2_pInt*o-1_pInt,pReal))/ &
-                       real(2_pInt*o-1_pInt,pReal)
-   dChemPot_dCh = dChemPot_dCh + 2.0_pReal*kBT*(2.0_pReal*Ch - 1.0_pReal)**real(2_pInt*o-2_pInt,pReal) 
- enddo
- 
- call hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCh, Ch, ip, el)
- ChemPot = ChemPot + KPot 
- dChemPot_dCh = dChemPot_dCh + dKPot_dCh
-
- if (Ch < 0.0_pReal) then
-   ChemPot = ChemPot - 3.0_pReal*hydrogenBoundPenalty*Ch*Ch
-   dChemPot_dCh = dChemPot_dCh - 6.0_pReal*hydrogenBoundPenalty*Ch
- elseif (Ch > 1.0_pReal) then
-   ChemPot = ChemPot + 3.0_pReal*hydrogenBoundPenalty*(1.0_pReal - Ch)*(1.0_pReal - Ch)
-   dChemPot_dCh = dChemPot_dCh - 6.0_pReal*hydrogenBoundPenalty*(1.0_pReal - Ch)
- endif       
- 
-end subroutine hydrogenflux_cahnhilliard_getChemPotAndItsTangent
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief updates hydrogen concentration with solution from Cahn-Hilliard PDE for solute transport
-!--------------------------------------------------------------------------------------------------
-subroutine hydrogenflux_cahnhilliard_putHydrogenConcAndItsRate(Ch,Chdot,ip,el)
- use material, only: &
-   mappingHomogenization, &
-   hydrogenConc, &
-   hydrogenConcRate, &
-   hydrogenfluxMapping
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   Ch, &
-   Chdot
- integer(pInt) :: &
-   homog, &
-   offset  
- 
- homog  = mappingHomogenization(2,ip,el)
- offset = hydrogenfluxMapping(homog)%p(ip,el)
- hydrogenConc    (homog)%p(offset) = Ch
- hydrogenConcRate(homog)%p(offset) = Chdot
-
-end subroutine hydrogenflux_cahnhilliard_putHydrogenConcAndItsRate
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of hydrogen transport results
-!--------------------------------------------------------------------------------------------------
-function hydrogenflux_cahnhilliard_postResults(ip,el)
- use material, only: &
-   mappingHomogenization, &
-   hydrogenflux_typeInstance, &
-   hydrogenConc, &
-   hydrogenfluxMapping
-
- implicit none
- integer(pInt),              intent(in) :: &
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal), dimension(hydrogenflux_cahnhilliard_sizePostResults(hydrogenflux_typeInstance(mappingHomogenization(2,ip,el)))) :: &
-   hydrogenflux_cahnhilliard_postResults
-
- integer(pInt) :: &
-   instance, homog, offset, o, c
-   
- homog     = mappingHomogenization(2,ip,el)
- offset    = hydrogenfluxMapping(homog)%p(ip,el)
- instance  = hydrogenflux_typeInstance(homog)
-
- c = 0_pInt
- hydrogenflux_cahnhilliard_postResults = 0.0_pReal
-
- do o = 1_pInt,hydrogenflux_cahnhilliard_Noutput(instance)
-    select case(hydrogenflux_cahnhilliard_outputID(o,instance))
- 
-      case (hydrogenConc_ID)
-        hydrogenflux_cahnhilliard_postResults(c+1_pInt) = hydrogenConc(homog)%p(offset)
-        c = c + 1
-    end select
- enddo
-end function hydrogenflux_cahnhilliard_postResults
-
-end module hydrogenflux_cahnhilliard
--- a/src/hydrogenflux_isoconc.f90
+++ b/src/hydrogenflux_isoconc.f90
@ -1,62 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for constant hydrogen concentration
-!--------------------------------------------------------------------------------------------------
-module hydrogenflux_isoconc
-
- implicit none
- private
- 
- public :: &
-   hydrogenflux_isoconc_init
-
-contains
-
-!--------------------------------------------------------------------------------------------------
-!> @brief allocates all neccessary fields, reads information from material configuration file
-!--------------------------------------------------------------------------------------------------
-subroutine hydrogenflux_isoconc_init()
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use prec, only: &
-   pReal, &
-   pInt 
- use IO, only: &
-   IO_timeStamp
- use material
- use config
- 
- implicit none
- integer(pInt) :: &
-   homog, &
-   NofMyHomog
-
- write(6,'(/,a)')   ' <<<+-  hydrogenflux_'//HYDROGENFLUX_isoconc_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
-
- initializeInstances: do homog = 1_pInt, material_Nhomogenization
-   
-   myhomog: if (hydrogenflux_type(homog) == HYDROGENFLUX_isoconc_ID) then
-     NofMyHomog = count(material_homog == homog)
-     hydrogenfluxState(homog)%sizeState = 0_pInt
-     hydrogenfluxState(homog)%sizePostResults = 0_pInt
-     allocate(hydrogenfluxState(homog)%state0   (0_pInt,NofMyHomog), source=0.0_pReal)
-     allocate(hydrogenfluxState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal)
-     allocate(hydrogenfluxState(homog)%state    (0_pInt,NofMyHomog), source=0.0_pReal)
-     
-     deallocate(hydrogenConc    (homog)%p)
-     deallocate(hydrogenConcRate(homog)%p)
-     allocate  (hydrogenConc    (homog)%p(1), source=hydrogenflux_initialCh(homog))
-     allocate  (hydrogenConcRate(homog)%p(1), source=0.0_pReal)
-
-   endif myhomog
- enddo initializeInstances
-
-
-end subroutine hydrogenflux_isoconc_init
-
-end module hydrogenflux_isoconc
--- a/src/kinematics_hydrogen_strain.f90
+++ b/src/kinematics_hydrogen_strain.f90
@ -1,263 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine incorporating kinematics resulting from interstitial hydrogen
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module kinematics_hydrogen_strain
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   kinematics_hydrogen_strain_sizePostResults, &                                                                !< cumulative size of post results
-   kinematics_hydrogen_strain_offset, &                                                                         !< which kinematics is my current damage mechanism?
-   kinematics_hydrogen_strain_instance                                                                          !< instance of damage kinematics mechanism
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
-   kinematics_hydrogen_strain_sizePostResult                                                                    !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
-   kinematics_hydrogen_strain_output                                                                            !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public  :: &
-   kinematics_hydrogen_strain_Noutput                                                                           !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),           allocatable,         private :: &
-   kinematics_hydrogen_strain_coeff
-
- public :: &
-   kinematics_hydrogen_strain_init, &
-   kinematics_hydrogen_strain_initialStrain, &
-   kinematics_hydrogen_strain_LiAndItsTangent, &
-   kinematics_hydrogen_strain_ChemPotAndItsTangent
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine kinematics_hydrogen_strain_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use debug, only: &
-   debug_level,&
-   debug_constitutive,&
-   debug_levelBasic
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   phase_kinematics, &
-   phase_Nkinematics, &
-   phase_Noutput, &
-   KINEMATICS_hydrogen_strain_label, &
-   KINEMATICS_hydrogen_strain_ID
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,phase,instance,kinematics
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
-
- maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt)
- if (maxNinstance == 0_pInt) return
- 
- if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
- 
- allocate(kinematics_hydrogen_strain_offset(material_Nphase), source=0_pInt)
- allocate(kinematics_hydrogen_strain_instance(material_Nphase), source=0_pInt)
- do phase = 1, material_Nphase
-   kinematics_hydrogen_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_hydrogen_strain_ID)
-   do kinematics = 1, phase_Nkinematics(phase)
-     if (phase_kinematics(kinematics,phase) == kinematics_hydrogen_strain_ID) &
-       kinematics_hydrogen_strain_offset(phase) = kinematics
-   enddo    
- enddo
-   
- allocate(kinematics_hydrogen_strain_sizePostResults(maxNinstance),                     source=0_pInt)
- allocate(kinematics_hydrogen_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
- allocate(kinematics_hydrogen_strain_output(maxval(phase_Noutput),maxNinstance))
-          kinematics_hydrogen_strain_output = ''
- allocate(kinematics_hydrogen_strain_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(kinematics_hydrogen_strain_coeff(maxNinstance),                               source=0.0_pReal) 
-
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_hydrogen_strain_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = kinematics_hydrogen_strain_instance(phase)                                                         ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('hydrogen_strain_coeff')
-         kinematics_hydrogen_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-     end select
-   endif; endif
- enddo parsingFile
-
-end subroutine kinematics_hydrogen_strain_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  report initial hydrogen strain based on current hydrogen conc deviation from 
-!>         equillibrium (0)
-!--------------------------------------------------------------------------------------------------
-pure function kinematics_hydrogen_strain_initialStrain(ipc, ip, el)
- use math, only: &
-   math_I3
- use material, only: &
-   material_phase, &
-   material_homog, &
-   hydrogenConc, &
-   hydrogenfluxMapping
- use lattice, only: &
-   lattice_equilibriumHydrogenConcentration
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   kinematics_hydrogen_strain_initialStrain                                                       !< initial thermal strain (should be small strain, though)
- integer(pInt) :: &
-   phase, &
-   homog, offset, instance
-   
- phase = material_phase(ipc,ip,el)
- instance = kinematics_hydrogen_strain_instance(phase)
- homog = material_homog(ip,el)
- offset = hydrogenfluxMapping(homog)%p(ip,el)
- 
- kinematics_hydrogen_strain_initialStrain = &
-   (hydrogenConc(homog)%p(offset) - lattice_equilibriumHydrogenConcentration(phase)) * &
-   kinematics_hydrogen_strain_coeff(instance)* math_I3
-  
-end function kinematics_hydrogen_strain_initialStrain
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  contains the constitutive equation for calculating the velocity gradient  
-!--------------------------------------------------------------------------------------------------
-subroutine kinematics_hydrogen_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
- use material, only: &
-   material_phase, &
-   material_homog, &
-   hydrogenConc, &
-   hydrogenConcRate, &
-   hydrogenfluxMapping
- use math, only: &
-   math_I3
- use lattice, only: &
-   lattice_equilibriumHydrogenConcentration
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(out), dimension(3,3) :: &
-   Li                                                                                               !< thermal velocity gradient
- real(pReal),   intent(out), dimension(3,3,3,3) :: &
-   dLi_dTstar3333                                                                                   !< derivative of Li with respect to Tstar (4th-order tensor)
- integer(pInt) :: &
-   phase, &
-   instance, &
-   homog, offset
- real(pReal) :: &
-   Ch, ChEq, ChDot  
-   
- phase = material_phase(ipc,ip,el)
- instance = kinematics_hydrogen_strain_instance(phase)
- homog = material_homog(ip,el)
- offset = hydrogenfluxMapping(homog)%p(ip,el)
- Ch = hydrogenConc(homog)%p(offset)
- ChDot = hydrogenConcRate(homog)%p(offset)
- ChEq = lattice_equilibriumHydrogenConcentration(phase)
- 
- Li = ChDot*math_I3* &
-      kinematics_hydrogen_strain_coeff(instance)/ &
-      (1.0_pReal + kinematics_hydrogen_strain_coeff(instance)*(Ch - ChEq))
- dLi_dTstar3333 = 0.0_pReal
-  
-end subroutine kinematics_hydrogen_strain_LiAndItsTangent
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  contains the kinematic contribution to hydrogen chemical potential  
-!--------------------------------------------------------------------------------------------------
-subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCh, Tstar_v, Fi0, Fi, ipc, ip, el)
- use material, only: &
-   material_phase
- use math, only: &
-   math_inv33, &
-   math_mul33x33, &
-   math_Mandel6to33, &
-   math_transpose33
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),  intent(in), dimension(6) :: &
-   Tstar_v
- real(pReal),  intent(in), dimension(3,3) :: &
-   Fi0, Fi
- real(pReal),   intent(out) :: &
-   ChemPot, dChemPot_dCh                                                               
- integer(pInt) :: &
-   phase, &
-   instance
-   
- phase = material_phase(ipc,ip,el)
- instance = kinematics_hydrogen_strain_instance(phase)
- 
- ChemPot = -kinematics_hydrogen_strain_coeff(instance)* &
-            sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* &
-                math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi))
- dChemPot_dCh = 0.0_pReal
-   
-end subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent
-
-end module kinematics_hydrogen_strain
--- a/src/kinematics_vacancy_strain.f90
+++ b/src/kinematics_vacancy_strain.f90
@ -1,264 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine incorporating kinematics resulting from vacancy point defects
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module kinematics_vacancy_strain
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   kinematics_vacancy_strain_sizePostResults, &                                                                !< cumulative size of post results
-   kinematics_vacancy_strain_offset, &                                                                         !< which kinematics is my current damage mechanism?
-   kinematics_vacancy_strain_instance                                                                          !< instance of damage kinematics mechanism
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
-   kinematics_vacancy_strain_sizePostResult                                                                    !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
-   kinematics_vacancy_strain_output                                                                            !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public  :: &
-   kinematics_vacancy_strain_Noutput                                                                           !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),           allocatable,         private :: &
-   kinematics_vacancy_strain_coeff
-
- public :: &
-   kinematics_vacancy_strain_init, &
-   kinematics_vacancy_strain_initialStrain, &
-   kinematics_vacancy_strain_LiAndItsTangent, &
-   kinematics_vacancy_strain_ChemPotAndItsTangent
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine kinematics_vacancy_strain_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use debug, only: &
-   debug_level,&
-   debug_constitutive,&
-   debug_levelBasic
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   phase_kinematics, &
-   phase_Nkinematics, &
-   phase_Noutput, &
-   KINEMATICS_vacancy_strain_label, &
-   KINEMATICS_vacancy_strain_ID
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,phase,instance,kinematics
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_vacancy_strain_LABEL//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
-
- maxNinstance = int(count(phase_kinematics == KINEMATICS_vacancy_strain_ID),pInt)
- if (maxNinstance == 0_pInt) return
- 
- if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
- 
- allocate(kinematics_vacancy_strain_offset(material_Nphase), source=0_pInt)
- allocate(kinematics_vacancy_strain_instance(material_Nphase), source=0_pInt)
- do phase = 1, material_Nphase
-   kinematics_vacancy_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_vacancy_strain_ID)
-   do kinematics = 1, phase_Nkinematics(phase)
-     if (phase_kinematics(kinematics,phase) == kinematics_vacancy_strain_ID) &
-       kinematics_vacancy_strain_offset(phase) = kinematics
-   enddo    
- enddo
-   
- allocate(kinematics_vacancy_strain_sizePostResults(maxNinstance),                     source=0_pInt)
- allocate(kinematics_vacancy_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
- allocate(kinematics_vacancy_strain_output(maxval(phase_Noutput),maxNinstance))
-          kinematics_vacancy_strain_output = ''
- allocate(kinematics_vacancy_strain_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(kinematics_vacancy_strain_coeff(maxNinstance),                               source=0.0_pReal) 
-
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_vacancy_strain_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = kinematics_vacancy_strain_instance(phase)                                                         ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('vacancy_strain_coeff')
-         kinematics_vacancy_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-     end select
-   endif; endif
- enddo parsingFile
-
-end subroutine kinematics_vacancy_strain_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  report initial vacancy strain based on current vacancy conc deviation from equillibrium
-!--------------------------------------------------------------------------------------------------
-pure function kinematics_vacancy_strain_initialStrain(ipc, ip, el)
- use math, only: &
-   math_I3
- use material, only: &
-   material_phase, &
-   material_homog, &
-   vacancyConc, &
-   vacancyfluxMapping
- use lattice, only: &
-   lattice_equilibriumVacancyConcentration
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   kinematics_vacancy_strain_initialStrain                                                       !< initial thermal strain (should be small strain, though)
- integer(pInt) :: &
-   phase, &
-   homog, offset, instance
-   
- phase = material_phase(ipc,ip,el)
- instance = kinematics_vacancy_strain_instance(phase)
- homog = material_homog(ip,el)
- offset = vacancyfluxMapping(homog)%p(ip,el)
- 
- kinematics_vacancy_strain_initialStrain = &
-   (vacancyConc(homog)%p(offset) - lattice_equilibriumVacancyConcentration(phase)) * &
-   kinematics_vacancy_strain_coeff(instance)* math_I3
-  
-end function kinematics_vacancy_strain_initialStrain
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  contains the constitutive equation for calculating the velocity gradient  
-!--------------------------------------------------------------------------------------------------
-subroutine kinematics_vacancy_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
- use material, only: &
-   material_phase, &
-   material_homog, &
-   vacancyConc, &
-   vacancyConcRate, &
-   vacancyfluxMapping
- use math, only: &
-   math_I3
- use lattice, only: &
-   lattice_equilibriumVacancyConcentration
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(out), dimension(3,3) :: &
-   Li                                                                                               !< thermal velocity gradient
- real(pReal),   intent(out), dimension(3,3,3,3) :: &
-   dLi_dTstar3333                                                                                   !< derivative of Li with respect to Tstar (4th-order tensor)
- integer(pInt) :: &
-   phase, &
-   instance, &
-   homog, offset
- real(pReal) :: &
-   Cv, CvEq, CvDot  
-   
- phase = material_phase(ipc,ip,el)
- instance = kinematics_vacancy_strain_instance(phase)
- homog = material_homog(ip,el)
- offset = vacancyfluxMapping(homog)%p(ip,el)
- 
- Cv = vacancyConc(homog)%p(offset)
- CvDot = vacancyConcRate(homog)%p(offset)
- CvEq = lattice_equilibriumvacancyConcentration(phase)
- 
- Li = CvDot*math_I3* &
-      kinematics_vacancy_strain_coeff(instance)/ &
-      (1.0_pReal + kinematics_vacancy_strain_coeff(instance)*(Cv - CvEq))
-      
- dLi_dTstar3333 = 0.0_pReal
-  
-end subroutine kinematics_vacancy_strain_LiAndItsTangent
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  contains the kinematic contribution to vacancy chemical potential  
-!--------------------------------------------------------------------------------------------------
-subroutine kinematics_vacancy_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCv, Tstar_v, Fi0, Fi, ipc, ip, el)
- use material, only: &
-   material_phase
- use math, only: &
-   math_inv33, &
-   math_mul33x33, &
-   math_Mandel6to33, &
-   math_transpose33
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),  intent(in), dimension(6) :: &
-   Tstar_v
- real(pReal),  intent(in), dimension(3,3) :: &
-   Fi0, Fi
- real(pReal),   intent(out) :: &
-   ChemPot, dChemPot_dCv                                                               
- integer(pInt) :: &
-   phase, &
-   instance
-   
- phase = material_phase(ipc,ip,el)
- instance = kinematics_vacancy_strain_instance(phase)
- 
- ChemPot = -kinematics_vacancy_strain_coeff(instance)* &
-            sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* &
-                math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi))
- dChemPot_dCv = 0.0_pReal
-   
-end subroutine kinematics_vacancy_strain_ChemPotAndItsTangent
-
-end module kinematics_vacancy_strain
--- a/src/material.f90
+++ b/src/material.f90
@ -36,29 +36,16 @@ module material
   SOURCE_damage_isoDuctile_label       = 'damage_isoductile', &
   SOURCE_damage_anisoBrittle_label     = 'damage_anisobrittle', &
   SOURCE_damage_anisoDuctile_label     = 'damage_anisoductile', &
-   SOURCE_vacancy_phenoplasticity_label = 'vacancy_phenoplasticity', &
-   SOURCE_vacancy_irradiation_label     = 'vacancy_irradiation', &
-   SOURCE_vacancy_thermalfluc_label     = 'vacancy_thermalfluctuation', &
   KINEMATICS_thermal_expansion_label   = 'thermal_expansion', &
   KINEMATICS_cleavage_opening_label    = 'cleavage_opening', &
   KINEMATICS_slipplane_opening_label   = 'slipplane_opening', &
-   KINEMATICS_vacancy_strain_label      = 'vacancy_strain', &
-   KINEMATICS_hydrogen_strain_label     = 'hydrogen_strain', &
   STIFFNESS_DEGRADATION_damage_label   = 'damage', &
-   STIFFNESS_DEGRADATION_porosity_label = 'porosity', &
   THERMAL_isothermal_label             = 'isothermal', &
   THERMAL_adiabatic_label              = 'adiabatic', &
   THERMAL_conduction_label             = 'conduction', &
   DAMAGE_none_label                    = 'none', &
   DAMAGE_local_label                   = 'local', &
   DAMAGE_nonlocal_label                = 'nonlocal', &
-   VACANCYFLUX_isoconc_label            = 'isoconcentration', &
-   VACANCYFLUX_isochempot_label         = 'isochemicalpotential', &
-   VACANCYFLUX_cahnhilliard_label       = 'cahnhilliard', &
-   POROSITY_none_label                  = 'none', &
-   POROSITY_phasefield_label            = 'phasefield', &
-   HYDROGENFLUX_isoconc_label           = 'isoconcentration', &
-   HYDROGENFLUX_cahnhilliard_label      = 'cahnhilliard', &
   HOMOGENIZATION_none_label            = 'none', &
   HOMOGENIZATION_isostrain_label       = 'isostrain', &
   HOMOGENIZATION_rgc_label             = 'rgc'
@ -87,25 +74,19 @@ module material
                 SOURCE_damage_isoBrittle_ID, &
                 SOURCE_damage_isoDuctile_ID, &
                 SOURCE_damage_anisoBrittle_ID, &
-                 SOURCE_damage_anisoDuctile_ID, &
-                 SOURCE_vacancy_phenoplasticity_ID, &
-                 SOURCE_vacancy_irradiation_ID, &
-                 SOURCE_vacancy_thermalfluc_ID
+                 SOURCE_damage_anisoDuctile_ID
 end enum

 enum, bind(c)
   enumerator :: KINEMATICS_undefined_ID, &
                 KINEMATICS_cleavage_opening_ID, &
                 KINEMATICS_slipplane_opening_ID, &
-                 KINEMATICS_thermal_expansion_ID, &
-                 KINEMATICS_vacancy_strain_ID, &
-                 KINEMATICS_hydrogen_strain_ID
+                 KINEMATICS_thermal_expansion_ID
 end enum

 enum, bind(c)
   enumerator :: STIFFNESS_DEGRADATION_undefined_ID, &
-                 STIFFNESS_DEGRADATION_damage_ID, &
-                 STIFFNESS_DEGRADATION_porosity_ID
+                 STIFFNESS_DEGRADATION_damage_ID
 end enum

 enum, bind(c)
@ -120,21 +101,6 @@ module material
                 DAMAGE_nonlocal_ID
 end enum

- enum, bind(c)
-   enumerator :: VACANCYFLUX_isoconc_ID, &
-                 VACANCYFLUX_isochempot_ID, &
-                 VACANCYFLUX_cahnhilliard_ID
- end enum
-
- enum, bind(c)
-   enumerator :: POROSITY_none_ID, &
-                 POROSITY_phasefield_ID
- end enum
- enum, bind(c)
-   enumerator :: HYDROGENFLUX_isoconc_ID, &
-                 HYDROGENFLUX_cahnhilliard_ID
- end enum
-
 enum, bind(c)
   enumerator :: HOMOGENIZATION_undefined_ID, &
                 HOMOGENIZATION_none_ID, &
@ -150,12 +116,6 @@ module material
   thermal_type                                                                                     !< thermal transport model
 integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable, public, protected :: &
   damage_type                                                                                      !< nonlocal damage model
- integer(kind(VACANCYFLUX_isoconc_ID)),      dimension(:),   allocatable, public, protected :: &
-   vacancyflux_type                                                                                 !< vacancy transport model
- integer(kind(POROSITY_none_ID)),            dimension(:),   allocatable, public, protected :: &
-   porosity_type                                                                                    !< porosity evolution model
- integer(kind(HYDROGENFLUX_isoconc_ID)),     dimension(:),   allocatable, public, protected :: &
-   hydrogenflux_type                                                                                !< hydrogen transport model

 integer(kind(SOURCE_undefined_ID)),         dimension(:,:), allocatable, public, protected :: &
   phase_source, &                                                                                  !< active sources mechanisms of each phase
@ -181,17 +141,11 @@ module material
   homogenization_typeInstance, &                                                                   !< instance of particular type of each homogenization
   thermal_typeInstance, &                                                                          !< instance of particular type of each thermal transport
   damage_typeInstance, &                                                                           !< instance of particular type of each nonlocal damage
-   vacancyflux_typeInstance, &                                                                      !< instance of particular type of each vacancy flux
-   porosity_typeInstance, &                                                                         !< instance of particular type of each porosity model
-   hydrogenflux_typeInstance, &                                                                     !< instance of particular type of each hydrogen flux
   microstructure_crystallite                                                                       !< crystallite setting ID of each microstructure ! DEPRECATED !!!!

 real(pReal), dimension(:), allocatable, public, protected :: &
   thermal_initialT, &                                                                              !< initial temperature per each homogenization
-   damage_initialPhi, &                                                                             !< initial damage per each homogenization
-   vacancyflux_initialCv, &                                                                         !< initial vacancy concentration per each homogenization
-   porosity_initialPhi, &                                                                           !< initial posority per each homogenization
-   hydrogenflux_initialCh                                                                           !< initial hydrogen concentration per each homogenization
+   damage_initialPhi                                                                                !< initial damage per each homogenization

 ! NEW MAPPINGS 
 integer(pInt), dimension(:), allocatable, public, protected :: &
@ -221,10 +175,7 @@ module material
 type(tState),        allocatable, dimension(:), public :: &
   homogState, &
   thermalState, &
-   damageState, &
-   vacancyfluxState, &
-   porosityState, &
-   hydrogenfluxState
+   damageState

 integer(pInt), dimension(:,:,:), allocatable, public, protected :: &
   material_texture                                                                                 !< texture (index) of each grain,IP,element
@ -274,20 +225,12 @@ module material

 type(tHomogMapping), allocatable, dimension(:), public :: &
   thermalMapping, &                                                                                !< mapping for thermal state/fields
-   damageMapping, &                                                                                 !< mapping for damage state/fields
-   vacancyfluxMapping, &                                                                            !< mapping for vacancy conc state/fields
-   porosityMapping, &                                                                               !< mapping for porosity state/fields
-   hydrogenfluxMapping                                                                              !< mapping for hydrogen conc state/fields
+   damageMapping                                                                                    !< mapping for damage state/fields

 type(group_float),  allocatable, dimension(:), public :: &
   temperature, &                                                                                   !< temperature field
   damage, &                                                                                        !< damage field
-   vacancyConc, &                                                                                   !< vacancy conc field
-   porosity, &                                                                                      !< porosity field
-   hydrogenConc, &                                                                                  !< hydrogen conc field
-   temperatureRate, &                                                                               !< temperature change rate field
-   vacancyConcRate, &                                                                               !< vacancy conc change field
-   hydrogenConcRate                                                                                 !< hydrogen conc change field
+   temperatureRate                                                                                  !< temperature change rate field

 public :: &
   material_init, &
@ -306,29 +249,16 @@ module material
   SOURCE_damage_isoDuctile_ID, &
   SOURCE_damage_anisoBrittle_ID, &
   SOURCE_damage_anisoDuctile_ID, &
-   SOURCE_vacancy_phenoplasticity_ID, &
-   SOURCE_vacancy_irradiation_ID, &
-   SOURCE_vacancy_thermalfluc_ID, &
   KINEMATICS_cleavage_opening_ID, &
   KINEMATICS_slipplane_opening_ID, &
   KINEMATICS_thermal_expansion_ID, &
-   KINEMATICS_vacancy_strain_ID, &
-   KINEMATICS_hydrogen_strain_ID, &
   STIFFNESS_DEGRADATION_damage_ID, &
-   STIFFNESS_DEGRADATION_porosity_ID, &
   THERMAL_isothermal_ID, &
   THERMAL_adiabatic_ID, &
   THERMAL_conduction_ID, &
   DAMAGE_none_ID, &
   DAMAGE_local_ID, &
   DAMAGE_nonlocal_ID, &
-   VACANCYFLUX_isoconc_ID, &
-   VACANCYFLUX_isochempot_ID, &
-   VACANCYFLUX_cahnhilliard_ID, &
-   POROSITY_none_ID, &
-   POROSITY_phasefield_ID, &
-   HYDROGENFLUX_isoconc_ID, &
-   HYDROGENFLUX_cahnhilliard_ID, &
   HOMOGENIZATION_none_ID, &
   HOMOGENIZATION_isostrain_ID, &
   HOMOGENIZATION_RGC_ID
@ -420,25 +350,14 @@ subroutine material_init()
 allocate(homogState         (size(config_homogenization)))
 allocate(thermalState       (size(config_homogenization)))
 allocate(damageState        (size(config_homogenization)))
- allocate(vacancyfluxState   (size(config_homogenization)))
- allocate(porosityState      (size(config_homogenization)))
- allocate(hydrogenfluxState  (size(config_homogenization)))

 allocate(thermalMapping     (size(config_homogenization)))
 allocate(damageMapping      (size(config_homogenization)))
- allocate(vacancyfluxMapping (size(config_homogenization)))
- allocate(porosityMapping    (size(config_homogenization)))
- allocate(hydrogenfluxMapping(size(config_homogenization)))

 allocate(temperature        (size(config_homogenization)))
 allocate(damage             (size(config_homogenization)))
- allocate(vacancyConc        (size(config_homogenization)))
- allocate(porosity           (size(config_homogenization)))
- allocate(hydrogenConc       (size(config_homogenization)))

 allocate(temperatureRate    (size(config_homogenization)))
- allocate(vacancyConcRate    (size(config_homogenization)))
- allocate(hydrogenConcRate   (size(config_homogenization)))

 do m = 1_pInt,size(config_microstructure)
   if(microstructure_crystallite(m) < 1_pInt .or. &
@ -511,17 +430,9 @@ subroutine material_init()
 do myHomog = 1,size(config_homogenization)
   thermalMapping     (myHomog)%p => mappingHomogenizationConst
   damageMapping      (myHomog)%p => mappingHomogenizationConst
-   vacancyfluxMapping (myHomog)%p => mappingHomogenizationConst
-   porosityMapping    (myHomog)%p => mappingHomogenizationConst
-   hydrogenfluxMapping(myHomog)%p => mappingHomogenizationConst
   allocate(temperature     (myHomog)%p(1), source=thermal_initialT(myHomog))
   allocate(damage          (myHomog)%p(1), source=damage_initialPhi(myHomog))
-   allocate(vacancyConc     (myHomog)%p(1), source=vacancyflux_initialCv(myHomog))
-   allocate(porosity        (myHomog)%p(1), source=porosity_initialPhi(myHomog))
-   allocate(hydrogenConc    (myHomog)%p(1), source=hydrogenflux_initialCh(myHomog))
   allocate(temperatureRate (myHomog)%p(1), source=0.0_pReal)
-   allocate(vacancyConcRate (myHomog)%p(1), source=0.0_pReal)
-   allocate(hydrogenConcRate(myHomog)%p(1), source=0.0_pReal)
 enddo

 end subroutine material_init
@ -545,23 +456,14 @@ subroutine material_parseHomogenization
 allocate(homogenization_type(size(config_homogenization)),           source=HOMOGENIZATION_undefined_ID)
 allocate(thermal_type(size(config_homogenization)),                  source=THERMAL_isothermal_ID)
 allocate(damage_type (size(config_homogenization)),                  source=DAMAGE_none_ID)
- allocate(vacancyflux_type(size(config_homogenization)),              source=VACANCYFLUX_isoconc_ID)
- allocate(porosity_type (size(config_homogenization)),                source=POROSITY_none_ID)
- allocate(hydrogenflux_type(size(config_homogenization)),             source=HYDROGENFLUX_isoconc_ID)
 allocate(homogenization_typeInstance(size(config_homogenization)),   source=0_pInt)
 allocate(thermal_typeInstance(size(config_homogenization)),          source=0_pInt)
 allocate(damage_typeInstance(size(config_homogenization)),           source=0_pInt)
- allocate(vacancyflux_typeInstance(size(config_homogenization)),      source=0_pInt)
- allocate(porosity_typeInstance(size(config_homogenization)),         source=0_pInt)
- allocate(hydrogenflux_typeInstance(size(config_homogenization)),     source=0_pInt)
 allocate(homogenization_Ngrains(size(config_homogenization)),        source=0_pInt)
 allocate(homogenization_Noutput(size(config_homogenization)),        source=0_pInt)
 allocate(homogenization_active(size(config_homogenization)),         source=.false.)  !!!!!!!!!!!!!!!
 allocate(thermal_initialT(size(config_homogenization)),              source=300.0_pReal)
 allocate(damage_initialPhi(size(config_homogenization)),             source=1.0_pReal)
- allocate(vacancyflux_initialCv(size(config_homogenization)),         source=0.0_pReal)
- allocate(porosity_initialPhi(size(config_homogenization)),           source=1.0_pReal)
- allocate(hydrogenflux_initialCh(size(config_homogenization)),        source=0.0_pReal)

 forall (h = 1_pInt:size(config_homogenization)) &
   homogenization_active(h) = any(mesh_homogenizationAt == h)
@ -620,53 +522,6 @@ subroutine material_parseHomogenization
     end select

   endif
-   
-   if (config_homogenization(h)%keyExists('vacancyflux')) then
-     vacancyflux_initialCv(h) = config_homogenization(h)%getFloat('cv0',defaultVal=0.0_pReal)
-
-     tag = config_homogenization(h)%getString('vacancyflux')
-     select case (trim(tag))
-       case(VACANCYFLUX_isoconc_label)
-         vacancyflux_type(h) = VACANCYFLUX_isoconc_ID
-       case(VACANCYFLUX_isochempot_label)
-         vacancyflux_type(h) = VACANCYFLUX_isochempot_ID
-       case(VACANCYFLUX_cahnhilliard_label)
-         vacancyflux_type(h) = VACANCYFLUX_cahnhilliard_ID
-       case default
-         call IO_error(500_pInt,ext_msg=trim(tag))
-     end select
-   
-   endif
-
-   if (config_homogenization(h)%keyExists('porosity')) then
-     !ToDo?
-
-     tag = config_homogenization(h)%getString('porosity')
-     select case (trim(tag))
-       case(POROSITY_NONE_label)
-         porosity_type(h) = POROSITY_none_ID
-       case(POROSITY_phasefield_label)
-         porosity_type(h) = POROSITY_phasefield_ID
-       case default
-         call IO_error(500_pInt,ext_msg=trim(tag))
-      end select
-
-   endif
-
-   if (config_homogenization(h)%keyExists('hydrogenflux')) then
-     hydrogenflux_initialCh(h) = config_homogenization(h)%getFloat('ch0',defaultVal=0.0_pReal)
-
-     tag = config_homogenization(h)%getString('hydrogenflux')
-     select case (trim(tag))
-       case(HYDROGENFLUX_isoconc_label)
-         hydrogenflux_type(h) = HYDROGENFLUX_isoconc_ID
-       case(HYDROGENFLUX_cahnhilliard_label)
-         hydrogenflux_type(h) = HYDROGENFLUX_cahnhilliard_ID
-       case default
-         call IO_error(500_pInt,ext_msg=trim(tag))
-     end select
-
-   endif

 enddo

@ -674,9 +529,6 @@ subroutine material_parseHomogenization
   homogenization_typeInstance(h)  = count(homogenization_type(1:h)  == homogenization_type(h))
   thermal_typeInstance(h)         = count(thermal_type       (1:h)  == thermal_type       (h))
   damage_typeInstance(h)          = count(damage_type        (1:h)  == damage_type        (h))
-   vacancyflux_typeInstance(h)     = count(vacancyflux_type   (1:h)  == vacancyflux_type   (h))
-   porosity_typeInstance(h)        = count(porosity_type      (1:h)  == porosity_type      (h))
-   hydrogenflux_typeInstance(h)    = count(hydrogenflux_type  (1:h)  == hydrogenflux_type  (h))
 enddo

 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active)
@ -866,12 +718,6 @@ subroutine material_parsePhase
         phase_source(sourceCtr,p) = SOURCE_damage_anisoBrittle_ID
       case (SOURCE_damage_anisoDuctile_label)
         phase_source(sourceCtr,p) = SOURCE_damage_anisoDuctile_ID
-       case (SOURCE_vacancy_phenoplasticity_label)
-         phase_source(sourceCtr,p) = SOURCE_vacancy_phenoplasticity_ID
-       case (SOURCE_vacancy_irradiation_label)
-         phase_source(sourceCtr,p) = SOURCE_vacancy_irradiation_ID
-       case (SOURCE_vacancy_thermalfluc_label)
-         phase_source(sourceCtr,p) = SOURCE_vacancy_thermalfluc_ID
     end select
   enddo

@ -890,10 +736,6 @@ subroutine material_parsePhase
         phase_kinematics(kinematicsCtr,p) = KINEMATICS_slipplane_opening_ID
       case (KINEMATICS_thermal_expansion_label)
         phase_kinematics(kinematicsCtr,p) = KINEMATICS_thermal_expansion_ID
-       case (KINEMATICS_vacancy_strain_label)
-         phase_kinematics(kinematicsCtr,p) = KINEMATICS_vacancy_strain_ID
-       case (KINEMATICS_hydrogen_strain_label)
-         phase_kinematics(kinematicsCtr,p) = KINEMATICS_hydrogen_strain_ID
     end select
   enddo
 #if defined(__GFORTRAN__)
@ -907,8 +749,6 @@ subroutine material_parsePhase
     select case (trim(str(stiffDegradationCtr)))
       case (STIFFNESS_DEGRADATION_damage_label)
         phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID
-       case (STIFFNESS_DEGRADATION_porosity_label)
-         phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_porosity_ID
    end select
   enddo
 enddo
--- a/src/plastic_isotropic.f90
+++ b/src/plastic_isotropic.f90
@ -1,8 +1,9 @@
 !--------------------------------------------------------------------------------------------------
 !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for isotropic (ISOTROPIC) plasticity
-!> @details Isotropic (ISOTROPIC) Plasticity which resembles the phenopowerlaw plasticity without
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for isotropic plasticity
+!> @details Isotropic Plasticity which resembles the phenopowerlaw plasticity without
 !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
 !! untextured polycrystal
 !--------------------------------------------------------------------------------------------------
@ -10,57 +11,58 @@ module plastic_isotropic
 use prec, only: &
   pReal,&
   pInt
- 
+
 implicit none
 private
 integer(pInt),                       dimension(:,:),   allocatable, target, public :: &
-   plastic_isotropic_sizePostResult                                                                   !< size of each post result output
+   plastic_isotropic_sizePostResult                                                                 !< size of each post result output
 character(len=64),                   dimension(:,:),   allocatable, target, public :: &
-   plastic_isotropic_output                                                                           !< name of each post result output
- integer(pInt),                       dimension(:),     allocatable, target, public :: &
-   plastic_isotropic_Noutput                                                                          !< number of outputs per instance
- 
- enum, bind(c) 
-   enumerator :: undefined_ID, &
-                 flowstress_ID, &
-                 strainrate_ID
+   plastic_isotropic_output                                                                         !< name of each post result output
+
+ enum, bind(c)
+   enumerator :: &
+     undefined_ID, &
+     flowstress_ID, &
+     strainrate_ID
 end enum

- type, private :: tParameters                                                                         !< container type for internal constitutive parameters
-   integer(kind(undefined_ID)), allocatable, dimension(:) :: & 
-     outputID
-  real(pReal) :: &
-     fTaylor, &
-     tau0, &
-     gdot0, &
-     n, &
+ type, private :: tParameters
+   real(pReal) :: &
+     fTaylor, &                                                                                     !< Taylor factor
+     tau0, &                                                                                        !< initial critical stress
+     gdot0, &                                                                                       !< reference strain rate
+     n, &                                                                                           !< stress exponent
     h0, &
     h0_slopeLnRate, &
-     tausat, &
+     tausat, &                                                                                      !< maximum critical stress
     a, &
-     aTolFlowstress, &
-     aTolShear, &
     tausat_SinhFitA, &
     tausat_SinhFitB, &
     tausat_SinhFitC, &
-     tausat_SinhFitD
-  logical :: &
+     tausat_SinhFitD, &
+     aTolFlowstress, &
+     aTolShear
+   integer(pInt) :: &
+     of_debug = 0_pInt
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: & 
+     outputID
+   logical :: &
     dilatation
 end type

- type(tParameters), dimension(:), allocatable, target, private :: param                               !< containers of constitutive parameters (len Ninstance)
- 
- type, private :: tIsotropicState                                                                     !< internal state aliases
-   real(pReal), pointer,     dimension(:) :: &                                                        ! scalars along NipcMyInstance
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
+
+ type, private :: tIsotropicState
+   real(pReal), pointer, dimension(:) :: &
     flowstress, &
     accumulatedShear
 end type

- type(tIsotropicState), allocatable, dimension(:), private :: &                                       !< state aliases per instance
-   state, &
-   dotState
+ type(tIsotropicState), allocatable, dimension(:), private :: &
+   dotState, &
+   state

- public  :: &
+ public :: &
   plastic_isotropic_init, &
   plastic_isotropic_LpAndItsTangent, &
   plastic_isotropic_LiAndItsTangent, &
@ -80,20 +82,29 @@ subroutine plastic_isotropic_init()
   compiler_version, &
   compiler_options
 #endif
-use IO
+ use prec, only: &
+   pStringLen
 use debug, only: &
+#ifdef DEBUG
+   debug_e, &
+   debug_i, &
+   debug_g, &
+   debug_levelExtensive, &
+#endif
   debug_level, &
-   debug_constitutive, &
+   debug_constitutive,&
   debug_levelBasic
- use numerics, only: &
-   numerics_integrator
- use math, only: &
-   math_Mandel3333to66, &
-   math_Voigt66to3333
+ use IO, only: &
+   IO_error, &
+   IO_timeStamp
 use material, only: &
+#ifdef DEBUG
+   phasememberAt, &
+#endif
   phase_plasticity, &
   phase_plasticityInstance, &
   phase_Noutput, &
+   material_allocatePlasticState, &
   PLASTICITY_ISOTROPIC_label, &
   PLASTICITY_ISOTROPIC_ID, &
   material_phase, &
@ -101,148 +112,140 @@ use IO
 use config, only: &
   MATERIAL_partPhase, &
   config_phase
-   
- use lattice  
+ use lattice

 implicit none
- 
- type(tParameters), pointer :: prm
- 
 integer(pInt) :: &
-   phase, & 
-   instance, &
-   maxNinstance, &
-   sizeDotState, &
-   sizeState, &
-   sizeDeltaState
- character(len=65536) :: &
-   extmsg    = ''
- integer(pInt) :: NipcMyPhase,i
- character(len=65536), dimension(:), allocatable :: outputs
+   Ninstance, &
+   p, i, &
+   NipcMyPhase, &
+   sizeState, sizeDotState
+   
+ character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
+ 
+ integer(kind(undefined_ID)) :: &
+   outputID

- write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'
+ character(len=pStringLen) :: &
+   extmsg = ''
+ character(len=65536), dimension(:), allocatable :: &
+   outputs
+
+ write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'
+ write(6,'(/,a)')   ' Maiti and Eisenlohr, Scripta Materialia, 145:37-40, 2018'
+ write(6,'(/,a)')   ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
 write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- 
- maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt)
+
+ Ninstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt)
 if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance

-! public variables
- allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
- allocate(plastic_isotropic_output(maxval(phase_Noutput), maxNinstance))
+ allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), Ninstance),source=0_pInt)
+ allocate(plastic_isotropic_output(maxval(phase_Noutput), Ninstance))
          plastic_isotropic_output = ''
- allocate(plastic_isotropic_Noutput(maxNinstance),                              source=0_pInt)

- allocate(param(maxNinstance))                                                                      ! one container of parameters per instance
- allocate(state(maxNinstance))                                                                      ! internal state aliases
- allocate(dotState(maxNinstance))
+ allocate(param(Ninstance))
+ allocate(state(Ninstance))
+ allocate(dotState(Ninstance))

- do phase = 1_pInt, size(phase_plasticityInstance)
-   if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then
-     instance = phase_plasticityInstance(phase)
-     prm => param(instance)                                                                         ! shorthand pointer to parameter object of my constitutive law
-     prm%tau0            =  config_phase(phase)%getFloat('tau0')
-     prm%tausat          =  config_phase(phase)%getFloat('tausat')
-     prm%gdot0           =  config_phase(phase)%getFloat('gdot0')
-     prm%n               =  config_phase(phase)%getFloat('n')
-     prm%h0              =  config_phase(phase)%getFloat('h0')
-     prm%fTaylor         =  config_phase(phase)%getFloat('m')
-     prm%h0_slopeLnRate  =  config_phase(phase)%getFloat('h0_slopelnrate', defaultVal=0.0_pReal)
-     prm%tausat_SinhFitA =  config_phase(phase)%getFloat('tausat_sinhfita',defaultVal=0.0_pReal)
-     prm%tausat_SinhFitB =  config_phase(phase)%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal)
-     prm%tausat_SinhFitC =  config_phase(phase)%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal)
-     prm%tausat_SinhFitD =  config_phase(phase)%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal)
-     prm%a               =  config_phase(phase)%getFloat('a')
-     prm%aTolFlowStress  =  config_phase(phase)%getFloat('atol_flowstress',defaultVal=1.0_pReal)
-     prm%aTolShear       =  config_phase(phase)%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
-    
-     prm%dilatation      = config_phase(phase)%keyExists('/dilatation/')
-
-#if defined(__GFORTRAN__)
-     outputs = ['GfortranBug86277']
-     outputs = config_phase(phase)%getStrings('(output)',defaultVal=outputs)
-     if (outputs(1) == 'GfortranBug86277') outputs = [character(len=65536)::]
-#else
-     outputs = config_phase(phase)%getStrings('(output)',defaultVal=[character(len=65536)::])
+ do p = 1_pInt, size(phase_plasticityInstance)
+   if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle
+   associate(prm => param(phase_plasticityInstance(p)), &
+             dot => dotState(phase_plasticityInstance(p)), &
+             stt => state(phase_plasticityInstance(p)), &
+             config => config_phase(p))
+   
+#ifdef DEBUG
+   if  (p==material_phase(debug_g,debug_i,debug_e)) then
+     prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
+   endif
 #endif
-     allocate(prm%outputID(0))
-     do i=1_pInt, size(outputs)
-       select case(outputs(i))
-         case ('flowstress')
-           plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
-           plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputs(i)
-           plasticState(phase)%sizePostResults =  plasticState(phase)%sizePostResults + 1_pInt
-           plastic_isotropic_sizePostResult(i,instance) = 1_pInt
-           prm%outputID = [prm%outputID,flowstress_ID]
-         case ('strainrate')
-           plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
-           plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputs(i)
-           plasticState(phase)%sizePostResults = &
-             plasticState(phase)%sizePostResults + 1_pInt
-           plastic_isotropic_sizePostResult(i,instance) = 1_pInt
-           prm%outputID = [prm%outputID,strainrate_ID]
-       end select
-     enddo
+
+   prm%tau0            =  config%getFloat('tau0')
+   prm%tausat          =  config%getFloat('tausat')
+   prm%gdot0           =  config%getFloat('gdot0')
+   prm%n               =  config%getFloat('n')
+   prm%h0              =  config%getFloat('h0')
+   prm%fTaylor         =  config%getFloat('m')
+   prm%h0_slopeLnRate  =  config%getFloat('h0_slopelnrate', defaultVal=0.0_pReal)
+   prm%tausat_SinhFitA =  config%getFloat('tausat_sinhfita',defaultVal=0.0_pReal)
+   prm%tausat_SinhFitB =  config%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal)
+   prm%tausat_SinhFitC =  config%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal)
+   prm%tausat_SinhFitD =  config%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal)
+   prm%a               =  config%getFloat('a')
+   prm%aTolFlowStress  =  config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
+   prm%aTolShear       =  config%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
+   
+   prm%dilatation      = config%keyExists('/dilatation/')

 !--------------------------------------------------------------------------------------------------
 !  sanity checks
-     extmsg = ''
-     if (prm%aTolShear        <= 0.0_pReal) extmsg = trim(extmsg)//"'aTolShear' "
-     if (prm%tau0              < 0.0_pReal) extmsg = trim(extmsg)//"'tau0' "
-     if (prm%gdot0            <= 0.0_pReal) extmsg = trim(extmsg)//"'gdot0' "
-     if (prm%n                <= 0.0_pReal) extmsg = trim(extmsg)//"'n' "
-     if (prm%tausat           <= prm%tau0)    extmsg = trim(extmsg)//"'tausat' "
-     if (prm%a                <= 0.0_pReal) extmsg = trim(extmsg)//"'a' " 
-     if (prm%fTaylor          <= 0.0_pReal) extmsg = trim(extmsg)//"'m' "
-     if (prm%aTolFlowstress   <= 0.0_pReal) extmsg = trim(extmsg)//"'atol_flowstress' "
-     if (extmsg /= '') call IO_error(211_pInt,ip=instance,&
-                            ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
+   extmsg = ''
+   if (prm%aTolShear        <= 0.0_pReal) extmsg = trim(extmsg)//'aTolShear ' 
+   if (prm%tau0              < 0.0_pReal) extmsg = trim(extmsg)//'tau0 ' 
+   if (prm%gdot0            <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0 ' 
+   if (prm%n                <= 0.0_pReal) extmsg = trim(extmsg)//'n ' 
+   if (prm%tausat           <= prm%tau0)  extmsg = trim(extmsg)//'tausat '
+   if (prm%a                <= 0.0_pReal) extmsg = trim(extmsg)//'a ' 
+   if (prm%fTaylor          <= 0.0_pReal) extmsg = trim(extmsg)//'m ' 
+   if (prm%aTolFlowstress   <= 0.0_pReal) extmsg = trim(extmsg)//'atol_flowstress '
+   if (prm%aTolShear        <= 0.0_pReal) extmsg = trim(extmsg)//'atol_shear ' 
+   
+!--------------------------------------------------------------------------------------------------
+!  exit if any parameter is out of range
+   if (extmsg /= '') &
+     call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
+
+!--------------------------------------------------------------------------------------------------
+!  output pararameters
+   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+   allocate(prm%outputID(0))
+   do i=1_pInt, size(outputs)
+     outputID = undefined_ID
+     select case(outputs(i))
+       case ('flowstress')
+         outputID = flowstress_ID
+       case ('strainrate')
+         outputID = strainrate_ID
+     end select
+     
+     if (outputID /= undefined_ID) then
+      plastic_isotropic_output(i,phase_plasticityInstance(p)) = outputs(i)
+      plastic_isotropic_sizePostResult(i,phase_plasticityInstance(p)) = 1_pInt
+      prm%outputID = [prm%outputID , outputID]
+    endif
+
+   end do

 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-     NipcMyPhase = count(material_phase == phase)                                                   ! number of own material points (including point components ipc)
+   NipcMyPhase = count(material_phase == p)
+   sizeState   = size(["flowstress       ","accumulated_shear"])
+   sizeDotState = sizeState

-     sizeDotState   = size(["flowstress       ","accumulated_shear"])
-     sizeDeltaState = 0_pInt                                                                         ! no sudden jumps in state
-     sizeState      = sizeDotState + sizeDeltaState
-     plasticState(phase)%sizeState = sizeState
-     plasticState(phase)%sizeDotState = sizeDotState
-     plasticState(phase)%sizeDeltaState = sizeDeltaState
-     plasticState(phase)%nSlip = 1
-     allocate(plasticState(phase)%aTolState          (   sizeState))
-     allocate(plasticState(phase)%state0             (   sizeState,NipcMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%partionedState0    (   sizeState,NipcMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%subState0          (   sizeState,NipcMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%state              (   sizeState,NipcMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%dotState           (sizeDotState,NipcMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState       (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
-       allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(plasticState(phase)%RK4dotState      (sizeDotState,NipcMyPhase),source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
+   call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, &
+                                      1_pInt,0_pInt,0_pInt)
+   plasticState(p)%sizePostResults = sum(plastic_isotropic_sizePostResult(:,phase_plasticityInstance(p)))

 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and aTolState
+   stt%flowstress  => plasticState(p)%state    (1,1:NipcMyPhase)
+   stt%flowstress  = prm%tau0
+   dot%flowstress  => plasticState(p)%dotState (1,1:NipcMyPhase)
+   plasticState(p)%aTolState(1)       =  prm%aTolFlowstress

-     state(instance)%flowstress             => plasticState(phase)%state    (1,1:NipcMyPhase)
-     dotState(instance)%flowstress          => plasticState(phase)%dotState (1,1:NipcMyPhase)
-     plasticState(phase)%state0(1,1:NipcMyPhase) = prm%tau0
-     plasticState(phase)%aTolState(1)       =  prm%aTolFlowstress
+   stt%accumulatedShear  => plasticState(p)%state    (2,1:NipcMyPhase)
+   dot%accumulatedShear  => plasticState(p)%dotState (2,1:NipcMyPhase)
+   plasticState(p)%aTolState(2)       =  prm%aTolShear
+   ! global alias
+   plasticState(p)%slipRate           => plasticState(p)%dotState(2:2,1:NipcMyPhase)
+   plasticState(p)%accumulatedSlip    => plasticState(p)%state   (2:2,1:NipcMyPhase)
+   
+   plasticState(p)%state0 = plasticState(p)%state                                                   ! ToDo: this could be done centrally
+   
+  end associate

-     state(instance)%accumulatedShear       => plasticState(phase)%state    (2,1:NipcMyPhase)
-     dotState(instance)%accumulatedShear    => plasticState(phase)%dotState (2,1:NipcMyPhase)
-     plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal
-     plasticState(phase)%aTolState(2)       =  prm%aTolShear
-     ! global alias
-     plasticState(phase)%slipRate           => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
-     plasticState(phase)%accumulatedSlip    => plasticState(phase)%state   (2:2,1:NipcMyPhase)
-
-endif
 enddo

 end subroutine plastic_isotropic_init
@ -251,308 +254,234 @@ end subroutine plastic_isotropic_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
-subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
+subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
+#ifdef DEBUG
 use debug, only: &
   debug_level, &
-   debug_constitutive, &
-   debug_levelBasic, &
+   debug_constitutive,&
   debug_levelExtensive, &
-   debug_levelSelective, &
-   debug_e, &
-   debug_i, &
-   debug_g
+   debug_levelSelective
+#endif
 use math, only: &
-   math_mul6x6, &
-   math_Mandel6to33, &
-   math_Plain3333to99, &
   math_deviatoric33, &
   math_mul33xx33
- use material, only: &
-   phasememberAt, &
-   material_phase, &
-   phase_plasticityInstance

 implicit none
- real(pReal), dimension(3,3), intent(out) :: &
+ real(pReal), dimension(3,3),     intent(out) :: &
   Lp                                                                                               !< plastic velocity gradient
- real(pReal), dimension(9,9), intent(out) :: &
-   dLp_dTstar99                                                                                     !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
+ real(pReal), dimension(3,3,3,3), intent(out) :: &
+   dLp_dMp                                                                                          !< derivative of Lp with respect to the Mandel stress

- real(pReal), dimension(6),   intent(in) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ real(pReal), dimension(3,3), intent(in) :: &
+   Mp                                                                                               !< Mandel stress
 integer(pInt),               intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
+   instance, &
+   of

- type(tParameters), pointer :: prm
- 
 real(pReal), dimension(3,3) :: &
-   Tstar_dev_33                                                                                     !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
- real(pReal), dimension(3,3,3,3) :: &
-   dLp_dTstar_3333                                                                                  !< derivative of Lp with respect to Tstar as 4th order tensor
+   Mp_dev                                                                                           !< deviatoric part of the Mandel stress
 real(pReal) :: &
   gamma_dot, &                                                                                     !< strainrate
-   norm_Tstar_dev, &                                                                                !< euclidean norm of Tstar_dev
-   squarenorm_Tstar_dev                                                                             !< square of the euclidean norm of Tstar_dev
+   norm_Mp_dev, &                                                                                   !< norm of the deviatoric part of the Mandel stress
+   squarenorm_Mp_dev                                                                                !< square of the norm of the deviatoric part of the Mandel stress
 integer(pInt) :: &
-   instance, of, &
   k, l, m, n

- of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
- prm => param(instance)
+ associate(prm => param(instance), stt => state(instance))
 
- Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v))                                        ! deviatoric part of 2nd Piola-Kirchhoff stress
- squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
- norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) 
+ Mp_dev = math_deviatoric33(Mp)
+ squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev)
+ norm_Mp_dev = sqrt(squarenorm_Mp_dev) 

- if (norm_Tstar_dev <= 0.0_pReal) then                                                              ! Tstar == 0 --> both Lp and dLp_dTstar are zero
-   Lp = 0.0_pReal
-   dLp_dTstar99 = 0.0_pReal
- else
-   gamma_dot = prm%gdot0 &
-             * ( sqrt(1.5_pReal) * norm_Tstar_dev / prm%fTaylor / state(instance)%flowstress(of) ) &
-             **prm%n
-
-   Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/prm%fTaylor 
+ if (norm_Mp_dev > 0.0_pReal) then
+   gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%fTaylor*stt%flowstress(of))) **prm%n

+   Lp = Mp_dev/norm_Mp_dev * gamma_dot/prm%fTaylor 
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
-       .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
-              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
-     write(6,'(a,i8,1x,i2,1x,i3)') '<< CONST isotropic >> at el ip g ',el,ip,ipc
+       .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
     write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Tstar (dev) / MPa', &
-                                      transpose(Tstar_dev_33(1:3,1:3))*1.0e-6_pReal
-     write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Tstar_dev*1.0e-6_pReal
+                                      transpose(Mp_dev)*1.0e-6_pReal
+     write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Mp_dev*1.0e-6_pReal
     write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', gamma_dot
   end if
-!--------------------------------------------------------------------------------------------------
-! Calculation of the tangent of Lp
+#endif
   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-     dLp_dTstar_3333(k,l,m,n) = (prm%n-1.0_pReal) * &
-                                      Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
+     dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
-     dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
+     dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal
   forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
-     dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
-   dLp_dTstar99 = math_Plain3333to99(gamma_dot / prm%fTaylor * &
-                                      dLp_dTstar_3333 / norm_Tstar_dev)
+     dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal
+   dLp_dMp = gamma_dot / prm%fTaylor * dLp_dMp / norm_Mp_dev
+ else
+   Lp = 0.0_pReal
+   dLp_dMp = 0.0_pReal
 end if
+ 
+ end associate
+ 
 end subroutine plastic_isotropic_LpAndItsTangent

+
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates plastic velocity gradient and its tangent
+! ToDo: Rename Tstar to Mi?
 !--------------------------------------------------------------------------------------------------
-subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,el)
+subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
 use math, only: &
-   math_mul6x6, &
-   math_Mandel6to33, &
-   math_Plain3333to99, &
   math_spherical33, &
   math_mul33xx33
- use material, only: &
-   phasememberAt, &
-   material_phase, &
-   phase_plasticityInstance

 implicit none
 real(pReal), dimension(3,3), intent(out) :: &
-   Li                                                                                               !< plastic velocity gradient
+   Li                                                                                               !< inleastic velocity gradient
 real(pReal), dimension(3,3,3,3), intent(out)  :: &
-   dLi_dTstar_3333                                                                                  !< derivative of Li with respect to Tstar as 4th order tensor
- real(pReal), dimension(6),   intent(in) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- integer(pInt),               intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
-
- type(tParameters), pointer :: prm
+   dLi_dTstar                                                                                       !< derivative of Li with respect to the Mandel stress
 
+ real(pReal), dimension(3,3),   intent(in) :: &
+   Tstar                                                                                            !< Mandel stress ToDo: Mi?
+ integer(pInt),                 intent(in) :: &
+   instance, &
+   of
+
 real(pReal), dimension(3,3) :: &
-   Tstar_sph_33                                                                                     !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
+   Tstar_sph                                                                                        !< sphiatoric part of the Mandel stress
 real(pReal) :: &
   gamma_dot, &                                                                                     !< strainrate
   norm_Tstar_sph, &                                                                                !< euclidean norm of Tstar_sph
   squarenorm_Tstar_sph                                                                             !< square of the euclidean norm of Tstar_sph
 integer(pInt) :: &
-   instance, of, &
   k, l, m, n

- of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
- prm => param(instance)
+ associate(prm => param(instance), stt => state(instance))
 
- Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v))                                         ! spherical part of 2nd Piola-Kirchhoff stress
- squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
+ Tstar_sph = math_spherical33(Tstar)
+ squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph,Tstar_sph)
 norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) 

- if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then                              ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
-   gamma_dot = prm%gdot0 &
-               * (sqrt(1.5_pReal) * norm_Tstar_sph / prm%fTaylor / state(instance)%flowstress(of) ) &
-               **prm%n
+ if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then                                          ! no stress or J2 plastitiy --> Li and its derivative are zero
+   gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%fTaylor*stt%flowstress(of))) **prm%n

-   Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/prm%fTaylor
-
-   !--------------------------------------------------------------------------------------------------
-   ! Calculation of the tangent of Li
+   Li = Tstar_sph/norm_Tstar_sph * gamma_dot/prm%fTaylor
   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-     dLi_dTstar_3333(k,l,m,n) = (prm%n-1.0_pReal) * &
-                                      Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
+     dLi_dTstar(k,l,m,n) = (prm%n-1.0_pReal) * Tstar_sph(k,l)*Tstar_sph(m,n) / squarenorm_Tstar_sph
   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
-     dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
+     dLi_dTstar(k,l,k,l) = dLi_dTstar(k,l,k,l) + 1.0_pReal

-   dLi_dTstar_3333 = gamma_dot / prm%fTaylor * &
-                                      dLi_dTstar_3333 / norm_Tstar_sph
+   dLi_dTstar = gamma_dot / prm%fTaylor * dLi_dTstar / norm_Tstar_sph
 else
-  Li = 0.0_pReal
-  dLi_dTstar_3333 = 0.0_pReal
+   Li = 0.0_pReal
+   dLi_dTstar = 0.0_pReal
 endif
+ 
+ end associate
+ 
 end subroutine plastic_isotropic_LiAndItsTangent


 !--------------------------------------------------------------------------------------------------
 !> @brief calculates the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
+subroutine plastic_isotropic_dotState(Mp,instance,of)
 use prec, only: &
   dEq0
 use math, only: &
-   math_mul6x6
- use material, only: &
-   phasememberAt, &
-   material_phase, &
-   phase_plasticityInstance
- 
+   math_mul33xx33, &
+   math_deviatoric33
+
 implicit none
- real(pReal), dimension(6), intent(in):: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- integer(pInt),             intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- type(tParameters), pointer :: prm
- real(pReal), dimension(6) :: &
-   Tstar_dev_v                                                                                      !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
+ real(pReal), dimension(3,3),  intent(in) :: &
+   Mp                                                                                               !< Mandel stress
+ integer(pInt),                intent(in) :: &
+   instance, &
+   of
+
 real(pReal) :: &
   gamma_dot, &                                                                                     !< strainrate
   hardening, &                                                                                     !< hardening coefficient
   saturation, &                                                                                    !< saturation flowstress
-   norm_Tstar_v                                                                                     !< euclidean norm of Tstar_dev
- integer(pInt) :: &
-   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
-   of                                                                                               !< shortcut notation for offset position in state array
+   norm_Mp                                                                                          !< norm of the (deviatoric) Mandel stress

- of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
- prm => param(instance)
+ associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
 
-!--------------------------------------------------------------------------------------------------
-! norm of (deviatoric) 2nd Piola-Kirchhoff stress
 if (prm%dilatation) then
-   norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
+   norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
 else
-   Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
-   Tstar_dev_v(4:6) = Tstar_v(4:6)
-   norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
- end if
-!--------------------------------------------------------------------------------------------------
-! strain rate 
- gamma_dot = prm%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & 
-            / &!-----------------------------------------------------------------------------------
-           (prm%fTaylor*state(instance)%flowstress(of) ))**prm%n
+   norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
+ endif
+ 
+ gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor*stt%flowstress(of))) **prm%n
 
-!--------------------------------------------------------------------------------------------------
-! hardening coefficient
 if (abs(gamma_dot) > 1e-12_pReal) then
   if (dEq0(prm%tausat_SinhFitA)) then
     saturation = prm%tausat
   else
     saturation = prm%tausat &
-                + asinh( (gamma_dot / prm%tausat_SinhFitA&
-                         )**(1.0_pReal / prm%tausat_SinhFitD)&
+                + asinh( (gamma_dot / prm%tausat_SinhFitA)**(1.0_pReal / prm%tausat_SinhFitD) &
                       )**(1.0_pReal / prm%tausat_SinhFitC) &
-                   / ( prm%tausat_SinhFitB &
-                       * (gamma_dot / prm%gdot0)**(1.0_pReal / prm%n) &
-                     )
+                   / prm%tausat_SinhFitB * (gamma_dot / prm%gdot0)**(1.0_pReal / prm%n)
   endif
   hardening = ( prm%h0 + prm%h0_slopeLnRate * log(gamma_dot) ) &
-               * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**prm%a &
-               * sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
+               * abs( 1.0_pReal - stt%flowstress(of)/saturation )**prm%a &
+               * sign(1.0_pReal, 1.0_pReal - stt%flowstress(of)/saturation)
 else
   hardening = 0.0_pReal
 endif

- dotState(instance)%flowstress      (of) = hardening * gamma_dot
- dotState(instance)%accumulatedShear(of) =             gamma_dot
+ dot%flowstress      (of) = hardening * gamma_dot
+ dot%accumulatedShear(of) =             gamma_dot
+ 
+ end associate

 end subroutine plastic_isotropic_dotState

+
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------
-function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
+function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
 use math, only: &
-   math_mul6x6
- use material, only: &
-   plasticState, &
-   material_phase, &
-   phasememberAt, &
-   phase_plasticityInstance
+   math_mul33xx33, &
+   math_deviatoric33

 implicit none
- real(pReal), dimension(6),  intent(in) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- integer(pInt),              intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
+ real(pReal), dimension(3,3),  intent(in) :: &
+   Mp                                                                                               !< Mandel stress
+ integer(pInt),                intent(in) :: &
+   instance, &
+   of

- type(tParameters), pointer :: prm
- 
- real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: &
-                                           plastic_isotropic_postResults
+ real(pReal), dimension(sum(plastic_isotropic_sizePostResult(:,instance))) :: &
+   postResults

- real(pReal), dimension(6) :: &
-   Tstar_dev_v                                                                                      !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
 real(pReal) :: &
-   norm_Tstar_v                                                                                     ! euclidean norm of Tstar_dev
+   norm_Mp                                                                                          !< norm of the Mandel stress
 integer(pInt) :: &
-   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
-   of, &                                                                                            !< shortcut notation for offset position in state array
-   c, &
-   o
+   o,c

- of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
- prm => param(instance)
+ associate(prm => param(instance), stt => state(instance))
 
-!--------------------------------------------------------------------------------------------------
-! norm of (deviatoric) 2nd Piola-Kirchhoff stress
 if (prm%dilatation) then
-   norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
+   norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
 else
-   Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
-   Tstar_dev_v(4:6) = Tstar_v(4:6)
-   norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
- end if
+   norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
+ endif
 
 c = 0_pInt
- plastic_isotropic_postResults = 0.0_pReal

- outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
+ outputsLoop: do o = 1_pInt,size(prm%outputID)
   select case(prm%outputID(o))
     case (flowstress_ID)
-       plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of)
+       postResults(c+1_pInt) = stt%flowstress(of)
       c = c + 1_pInt
     case (strainrate_ID)
-       plastic_isotropic_postResults(c+1_pInt) = &
-                prm%gdot0 * (            sqrt(1.5_pReal) * norm_Tstar_v & 
-             / &!----------------------------------------------------------------------------------
-              (prm%fTaylor * state(instance)%flowstress(of)) ) ** prm%n
+       postResults(c+1_pInt) = prm%gdot0 &
+                             * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor * stt%flowstress(of)))**prm%n
       c = c + 1_pInt
   end select
 enddo outputsLoop
+ 
+ end associate

 end function plastic_isotropic_postResults

--- a/src/plastic_none.f90
+++ b/src/plastic_none.f90
@ -48,7 +48,7 @@ subroutine plastic_none_init
   phase, &
   NofMyPhase
 
- write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_NONE_label//' init  -+>>>'
+ write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_NONE_label//' init  -+>>>'
 write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
 
--- a/src/plastic_phenopowerlaw.f90
+++ b/src/plastic_phenopowerlaw.f90
@ -82,8 +82,7 @@ module plastic_phenopowerlaw
     xi_slip, &
     xi_twin, &
     gamma_slip, &
-     gamma_twin, &
-     whole
+     gamma_twin
 end type

 type(tPhenopowerlawState), allocatable, dimension(:), private :: &
@ -95,6 +94,9 @@ module plastic_phenopowerlaw
   plastic_phenopowerlaw_LpAndItsTangent, &
   plastic_phenopowerlaw_dotState, &
   plastic_phenopowerlaw_postResults
+ private :: &
+   kinetics_slip, &
+   kinetics_twin

 contains

@ -110,8 +112,7 @@ subroutine plastic_phenopowerlaw_init
   compiler_options
 #endif
 use prec, only: &
-   pStringLen, &
-   dEq0
+   pStringLen
 use debug, only: &
   debug_level, &
   debug_constitutive,&
@ -119,7 +120,6 @@ subroutine plastic_phenopowerlaw_init
 use math, only: &
   math_expand
 use IO, only: &
-   IO_warning, &
   IO_error, &
   IO_timeStamp
 use material, only: &
@ -149,7 +149,7 @@ subroutine plastic_phenopowerlaw_init
 character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]

 integer(kind(undefined_ID)) :: &
-   outputID                                                                                         !< ID of each post result output
+   outputID

 character(len=pStringLen) :: &
   structure = '',&
@ -157,7 +157,7 @@ subroutine plastic_phenopowerlaw_init
 character(len=65536), dimension(:), allocatable :: &
   outputs

- write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init  -+>>>'
+ write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init  -+>>>'
 write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"

@ -177,20 +177,21 @@ subroutine plastic_phenopowerlaw_init
   if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
   associate(prm => param(phase_plasticityInstance(p)), &
             dot => dotState(phase_plasticityInstance(p)), &
-             stt => state(phase_plasticityInstance(p)))
+             stt => state(phase_plasticityInstance(p)), &
+             config => config_phase(p))

-   structure          = config_phase(p)%getString('lattice_structure')
+   structure          = config%getString('lattice_structure')

 !--------------------------------------------------------------------------------------------------
 !  optional parameters that need to be defined
-   prm%twinB          = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal)
-   prm%twinC          = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal)
-   prm%twinD          = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal)
-   prm%twinE          = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal)
+   prm%twinB          = config%getFloat('twin_b',defaultVal=1.0_pReal)
+   prm%twinC          = config%getFloat('twin_c',defaultVal=0.0_pReal)
+   prm%twinD          = config%getFloat('twin_d',defaultVal=0.0_pReal)
+   prm%twinE          = config%getFloat('twin_e',defaultVal=0.0_pReal)

-   prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal)
-   prm%aTolShear      = config_phase(p)%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
-   prm%aTolTwinfrac   = config_phase(p)%getFloat('atol_twinfrac',  defaultVal=1.0e-6_pReal)
+   prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
+   prm%aTolShear      = config%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
+   prm%aTolTwinfrac   = config%getFloat('atol_twinfrac',  defaultVal=1.0e-6_pReal)
   
   ! sanity checks
   if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//'aTolresistance '
@ -199,14 +200,14 @@ subroutine plastic_phenopowerlaw_init

 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
-   prm%Nslip      = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray)
+   prm%Nslip      = config%getInts('nslip',defaultVal=emptyIntArray)
   prm%totalNslip = sum(prm%Nslip)
   slipActive: if (prm%totalNslip > 0_pInt) then
     prm%Schmid_slip          = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
-                                                          config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
+                                                          config%getFloat('c/a',defaultVal=0.0_pReal))
     if(structure=='bcc') then
-       prm%nonSchmidCoeff     = config_phase(p)%getFloats('nonschmid_coefficients',&
-                                                          defaultVal = emptyRealArray)
+       prm%nonSchmidCoeff     = config%getFloats('nonschmid_coefficients',&
+                                                 defaultVal = emptyRealArray)
       prm%nonSchmid_pos      = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
       prm%nonSchmid_neg      = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
     else
@ -214,18 +215,18 @@ subroutine plastic_phenopowerlaw_init
       prm%nonSchmid_neg      = prm%Schmid_slip
     endif
     prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
-                                                             config_phase(p)%getFloats('interaction_slipslip'), &
+                                                             config%getFloats('interaction_slipslip'), &
                                                             structure(1:3))

-     prm%xi_slip_0            = config_phase(p)%getFloats('tau0_slip',   requiredShape=shape(prm%Nslip))
-     prm%xi_slip_sat          = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip))
-     prm%H_int                = config_phase(p)%getFloats('h_int',       requiredShape=shape(prm%Nslip), &
-                                                                         defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
+     prm%xi_slip_0            = config%getFloats('tau0_slip',   requiredSize=size(prm%Nslip))
+     prm%xi_slip_sat          = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
+     prm%H_int                = config%getFloats('h_int',       requiredSize=size(prm%Nslip), &
+                                                                defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
 
-     prm%gdot0_slip           = config_phase(p)%getFloat('gdot0_slip')
-     prm%n_slip               = config_phase(p)%getFloat('n_slip')
-     prm%a_slip               = config_phase(p)%getFloat('a_slip')
-     prm%h0_SlipSlip          = config_phase(p)%getFloat('h0_slipslip')
+     prm%gdot0_slip           = config%getFloat('gdot0_slip')
+     prm%n_slip               = config%getFloat('n_slip')
+     prm%a_slip               = config%getFloat('a_slip')
+     prm%h0_SlipSlip          = config%getFloat('h0_slipslip')

     ! expand: family => system
     prm%xi_slip_0   = math_expand(prm%xi_slip_0,  prm%Nslip)
@ -233,9 +234,9 @@ subroutine plastic_phenopowerlaw_init
     prm%H_int       = math_expand(prm%H_int,      prm%Nslip)

     ! sanity checks
-     if (prm%gdot0_slip <= 0.0_pReal)           extmsg = trim(extmsg)//'gdot0_slip '
-     if (dEq0(prm%a_slip))                      extmsg = trim(extmsg)//'a_slip '                    ! ToDo: negative values ok?
-     if (dEq0(prm%n_slip))                      extmsg = trim(extmsg)//'n_slip '                    ! ToDo: negative values ok?
+     if (prm%gdot0_slip      <= 0.0_pReal)      extmsg = trim(extmsg)//'gdot0_slip '
+     if (prm%a_slip          <= 0.0_pReal)      extmsg = trim(extmsg)//'a_slip '
+     if (prm%n_slip          <= 0.0_pReal)      extmsg = trim(extmsg)//'n_slip '
     if (any(prm%xi_slip_0   <= 0.0_pReal))     extmsg = trim(extmsg)//'xi_slip_0 '
     if (any(prm%xi_slip_sat <  prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat '
   else slipActive
@ -245,30 +246,30 @@ subroutine plastic_phenopowerlaw_init

 !--------------------------------------------------------------------------------------------------
 ! twin related parameters
-   prm%Ntwin      = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray)
+   prm%Ntwin      = config%getInts('ntwin', defaultVal=emptyIntArray)
   prm%totalNtwin = sum(prm%Ntwin)
   twinActive: if (prm%totalNtwin > 0_pInt) then
     prm%Schmid_twin          = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),&
-                                                          config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
+                                                          config%getFloat('c/a',defaultVal=0.0_pReal))
     prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
-                                                             config_phase(p)%getFloats('interaction_twintwin'), &
+                                                             config%getFloats('interaction_twintwin'), &
                                                             structure(1:3))
     prm%gamma_twin_char      = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),&
-                                                            config_phase(p)%getFloat('c/a'))
+                                                            config%getFloat('c/a'))

-     prm%xi_twin_0            = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin))
+     prm%xi_twin_0            = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))

-     prm%gdot0_twin           = config_phase(p)%getFloat('gdot0_twin')
-     prm%n_twin               = config_phase(p)%getFloat('n_twin')
-     prm%spr                  = config_phase(p)%getFloat('s_pr')
-     prm%h0_TwinTwin          = config_phase(p)%getFloat('h0_twintwin')
+     prm%gdot0_twin           = config%getFloat('gdot0_twin')
+     prm%n_twin               = config%getFloat('n_twin')
+     prm%spr                  = config%getFloat('s_pr')
+     prm%h0_TwinTwin          = config%getFloat('h0_twintwin')

     ! expand: family => system
     prm%xi_twin_0   = math_expand(prm%xi_twin_0,  prm%Ntwin)

     ! sanity checks
     if (prm%gdot0_twin <= 0.0_pReal)  extmsg = trim(extmsg)//'gdot0_twin '
-     if (dEq0(prm%n_twin))             extmsg = trim(extmsg)//'n_twin '                             ! ToDo: negative values ok?
+     if (prm%n_twin     <= 0.0_pReal)  extmsg = trim(extmsg)//'n_twin '
   else twinActive
     allocate(prm%interaction_TwinTwin(0,0))
     allocate(prm%xi_twin_0(0))
@ -279,10 +280,10 @@ subroutine plastic_phenopowerlaw_init
 ! slip-twin related parameters
   slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
     prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
-                                                             config_phase(p)%getFloats('interaction_sliptwin'), &
+                                                             config%getFloats('interaction_sliptwin'), &
                                                             structure(1:3))
     prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
-                                                             config_phase(p)%getFloats('interaction_twinslip'), &
+                                                             config%getFloats('interaction_twinslip'), &
                                                             structure(1:3))
   else slipAndTwinActive
     allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin))                              ! at least one dimension is 0
@ -297,7 +298,7 @@ subroutine plastic_phenopowerlaw_init

 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-   outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray)
+   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
   allocate(prm%outputID(0))
   do i=1_pInt, size(outputs)
     outputID = undefined_ID
@ -340,7 +341,7 @@ subroutine plastic_phenopowerlaw_init

 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-   NipcMyPhase = count(material_phase == p)                                                         ! number of IPCs containing my phase
+   NipcMyPhase = count(material_phase == p)
   sizeState = size(['tau_slip  ','gamma_slip']) * prm%TotalNslip &
             + size(['tau_twin  ','gamma_twin']) * prm%TotalNtwin
   sizeDotState = sizeState
@ -349,7 +350,6 @@ subroutine plastic_phenopowerlaw_init
                                      prm%totalNslip,prm%totalNtwin,0_pInt)
   plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,phase_plasticityInstance(p)))

-
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and aTolState
   startIndex = 1_pInt
@ -381,8 +381,7 @@ subroutine plastic_phenopowerlaw_init
   dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
   plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear

-   plasticState(p)%state0 = plasticState(p)%state                                                 ! ToDo: this could be done centrally
-   dot%whole => plasticState(p)%dotState
+   plasticState(p)%state0 = plasticState(p)%state                                                   ! ToDo: this could be done centrally

   end associate
 enddo
@ -398,7 +397,7 @@ end subroutine plastic_phenopowerlaw_init
 subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)

 implicit none
- real(pReal), dimension(3,3), intent(out) :: &
+ real(pReal), dimension(3,3),     intent(out) :: &
   Lp                                                                                               !< plastic velocity gradient
 real(pReal), dimension(3,3,3,3), intent(out) :: &
   dLp_dMp                                                                                          !< derivative of Lp with respect to the Mandel stress
@ -420,9 +419,9 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
 Lp = 0.0_pReal
 dLp_dMp = 0.0_pReal
 
- associate(prm => param(instance), stt => state(instance))
+ associate(prm => param(instance))
 
- call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
+ call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
 slipSystems: do i = 1_pInt, prm%totalNslip
   Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
@ -431,7 +430,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
                      + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
 enddo slipSystems

- call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin)
+ call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
 twinSystems: do i = 1_pInt, prm%totalNtwin
   Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
@ -452,7 +451,7 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
 implicit none
 real(pReal), dimension(3,3),  intent(in) :: &
   Mp                                                                                               !< Mandel stress
- integer(pInt),              intent(in) :: &
+ integer(pInt),                intent(in) :: &
   instance, &
   of

@ -468,7 +467,6 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)

 associate(prm => param(instance),  stt => state(instance),  dot => dotState(instance))

- dot%whole(:,of) = 0.0_pReal
 sumGamma = sum(stt%gamma_slip(:,of))
 sumF     = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)

@ -487,9 +485,9 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)

 !--------------------------------------------------------------------------------------------------
 ! shear rates
- call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg)
+ call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
 dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
- call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of))
+ call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of))

 !--------------------------------------------------------------------------------------------------
 ! hardening
@ -510,40 +508,108 @@ end subroutine plastic_phenopowerlaw_dotState


 !--------------------------------------------------------------------------------------------------
-!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
-!> @details Shear rates are calculated only optionally. 
+!> @brief return array of constitutive results
+!--------------------------------------------------------------------------------------------------
+function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
+ use math, only: &
+   math_mul33xx33
+
+ implicit none
+ real(pReal), dimension(3,3),  intent(in) :: &
+   Mp                                                                                               !< Mandel stress
+ integer(pInt),                intent(in) :: &
+   instance, &
+   of
+
+ real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: &
+   postResults
+
+ integer(pInt) :: &
+   o,c,i
+ real(pReal), dimension(param(instance)%totalNslip) :: &
+   gdot_slip_pos,gdot_slip_neg
+
+ c = 0_pInt
+ 
+ associate( prm => param(instance), stt => state(instance))
+
+ outputsLoop: do o = 1_pInt,size(prm%outputID)
+   select case(prm%outputID(o))
+
+     case (resistance_slip_ID)
+       postResults(c+1_pInt:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of)
+       c = c + prm%totalNslip
+     case (accumulatedshear_slip_ID)
+       postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of)
+       c = c + prm%totalNslip
+     case (shearrate_slip_ID)
+       call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
+       postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg
+       c = c + prm%totalNslip
+     case (resolvedstress_slip_ID)
+       do i = 1_pInt, prm%totalNslip
+         postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i))
+       enddo
+       c = c + prm%totalNslip
+
+     case (resistance_twin_ID)
+       postResults(c+1_pInt:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of)
+       c = c + prm%totalNtwin
+     case (accumulatedshear_twin_ID)
+       postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of)
+       c = c + prm%totalNtwin
+     case (shearrate_twin_ID)
+       call kinetics_twin(Mp,instance,of,postResults(c+1_pInt:c+prm%totalNtwin))
+       c = c + prm%totalNtwin
+     case (resolvedstress_twin_ID)
+       do i = 1_pInt, prm%totalNtwin
+         postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
+       enddo
+       c = c + prm%totalNtwin
+
+   end select
+ enddo outputsLoop
+ 
+ end associate
+
+end function plastic_phenopowerlaw_postResults
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
+!> @details Derivatives are calculated only optionally. 
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to 
 ! have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, &
-                           dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
+pure subroutine kinetics_slip(Mp,instance,of, &
+                              gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
 use prec, only: &
  dNeq0
 use math, only: &
   math_mul33xx33

 implicit none
- type(tParameters), intent(in) :: &
-   prm
- type(tPhenopowerlawState), intent(in) :: &
-   stt
- integer(pInt),     intent(in) :: &
+ real(pReal), dimension(3,3),  intent(in) :: &
+   Mp                                                                                               !< Mandel stress
+ integer(pInt),                intent(in) :: &
+   instance, &
   of
- real(pReal), dimension(prm%totalNslip), intent(out) :: &
+   
+ real(pReal), dimension(param(instance)%totalNslip), intent(out) :: &
   gdot_slip_pos, &
   gdot_slip_neg
- real(pReal), dimension(prm%totalNslip), optional, intent(out) :: &
+ real(pReal), dimension(param(instance)%totalNslip), intent(out), optional :: &
   dgdot_dtau_slip_pos, &
   dgdot_dtau_slip_neg
- real(pReal), dimension(3,3), intent(in) :: &
-   Mp

- real(pReal), dimension(prm%totalNslip) :: &
+ real(pReal), dimension(param(instance)%totalNslip) :: &
   tau_slip_pos, &
   tau_slip_neg
 integer(pInt) :: i
 logical       :: nonSchmidActive

+ associate(prm => param(instance), stt => state(instance))
+ 
 nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt

 do i = 1_pInt, prm%totalNslip
@ -580,40 +646,42 @@ pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, &
     dgdot_dtau_slip_neg = 0.0_pReal
   end where
 endif
+ end associate

 end subroutine kinetics_slip


 !--------------------------------------------------------------------------------------------------
-!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress.
+!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
 !  twinning is assumed to take place only in untwinned volume.
 !> @details Derivates are calculated only optionally.
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to 
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin)
+pure subroutine kinetics_twin(Mp,instance,of,&
+                              gdot_twin,dgdot_dtau_twin)
 use prec, only: &
  dNeq0
 use math, only: &
   math_mul33xx33

 implicit none
- type(tParameters), intent(in) :: &
-   prm
- type(tPhenopowerlawState), intent(in) :: &
-   stt
- integer(pInt),     intent(in) :: &
+ real(pReal), dimension(3,3),  intent(in) :: &
+   Mp                                                                                               !< Mandel stress
+ integer(pInt),                intent(in) :: &
+   instance, &
   of
- real(pReal), dimension(3,3), intent(in) :: &
-   Mp
- real(pReal), dimension(prm%totalNtwin), intent(out) :: &
+   
+ real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
   gdot_twin
- real(pReal), dimension(prm%totalNtwin), optional, intent(out) :: &
+ real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
   dgdot_dtau_twin

- real(pReal), dimension(prm%totalNtwin) :: &
+ real(pReal), dimension(param(instance)%totalNtwin) :: &
   tau_twin
 integer(pInt) :: i
+ 
+ associate(prm => param(instance), stt => state(instance))

 do i = 1_pInt, prm%totalNtwin
   tau_twin(i)  = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
@ -633,76 +701,9 @@ pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin)
     dgdot_dtau_twin = 0.0_pReal
   end where
 endif
-
-end subroutine kinetics_twin
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of constitutive results
-!--------------------------------------------------------------------------------------------------
-function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
- use math, only: &
-   math_mul33xx33
-
- implicit none
- real(pReal), dimension(3,3), intent(in) :: &
-   Mp                                                                                               !< Mandel stress
- integer(pInt),             intent(in) :: &
-   instance, &
-   of
-
- real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: &
-   postResults
-
- integer(pInt) :: &
-   o,c,i
- real(pReal), dimension(param(instance)%totalNslip) :: &
-   gdot_slip_pos,gdot_slip_neg
-
- postResults = 0.0_pReal
- c = 0_pInt
- 
- associate( prm => param(instance), stt => state(instance))
-
- outputsLoop: do o = 1_pInt,size(prm%outputID)
-   select case(prm%outputID(o))
-
-     case (resistance_slip_ID)
-       postResults(c+1_pInt:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of)
-       c = c + prm%totalNslip
-     case (accumulatedshear_slip_ID)
-       postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of)
-       c = c + prm%totalNslip
-     case (shearrate_slip_ID)
-       call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg)
-       postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg
-       c = c + prm%totalNslip
-     case (resolvedstress_slip_ID)
-       do i = 1_pInt, prm%totalNslip
-         postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i))
-       enddo
-       c = c + prm%totalNslip
-
-     case (resistance_twin_ID)
-       postResults(c+1_pInt:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of)
-       c = c + prm%totalNtwin
-     case (accumulatedshear_twin_ID)
-       postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of)
-       c = c + prm%totalNtwin
-     case (shearrate_twin_ID)
-       call kinetics_twin(prm,stt,of,Mp,postResults(c+1_pInt:c+prm%totalNtwin))
-       c = c + prm%totalNtwin
-     case (resolvedstress_twin_ID)
-       do i = 1_pInt, prm%totalNtwin
-         postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
-       enddo
-       c = c + prm%totalNtwin
-
-   end select
- enddo outputsLoop
 
 end associate

-end function plastic_phenopowerlaw_postResults
+end subroutine kinetics_twin

 end module plastic_phenopowerlaw
--- a/src/porosity_none.f90
+++ b/src/porosity_none.f90
@ -1,60 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for constant porosity
-!--------------------------------------------------------------------------------------------------
-module porosity_none
-
- implicit none
- private
- 
- public :: &
-   porosity_none_init
-
-contains
-
-!--------------------------------------------------------------------------------------------------
-!> @brief allocates all neccessary fields, reads information from material configuration file
-!--------------------------------------------------------------------------------------------------
-subroutine porosity_none_init()
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use prec, only: &
-   pReal, &
-   pInt 
- use IO, only: &
-   IO_timeStamp
- use material
- use config
- 
- implicit none
- integer(pInt) :: &
-   homog, &
-   NofMyHomog
-
- write(6,'(/,a)')   ' <<<+-  porosity_'//POROSITY_none_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
-
- initializeInstances: do homog = 1_pInt, material_Nhomogenization
-   
-   myhomog: if (porosity_type(homog) == POROSITY_none_ID) then
-     NofMyHomog = count(material_homog == homog)
-     porosityState(homog)%sizeState = 0_pInt
-     porosityState(homog)%sizePostResults = 0_pInt
-     allocate(porosityState(homog)%state0   (0_pInt,NofMyHomog), source=0.0_pReal)
-     allocate(porosityState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal)
-     allocate(porosityState(homog)%state    (0_pInt,NofMyHomog), source=0.0_pReal)
-     
-     deallocate(porosity(homog)%p)
-     allocate  (porosity(homog)%p(1), source=porosity_initialPhi(homog))
-     
-   endif myhomog
- enddo initializeInstances
-
-
-end subroutine porosity_none_init
-
-end module porosity_none
--- a/src/porosity_phasefield.f90
+++ b/src/porosity_phasefield.f90
@ -1,448 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for phase field modelling of pore nucleation and growth
-!> @details phase field model for pore nucleation and growth based on vacancy clustering
-!--------------------------------------------------------------------------------------------------
-module porosity_phasefield
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   porosity_phasefield_sizePostResults                                                           !< cumulative size of post results
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   porosity_phasefield_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   porosity_phasefield_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   porosity_phasefield_Noutput                                                                   !< number of outputs per instance of this porosity 
-
- enum, bind(c) 
-   enumerator :: undefined_ID, &
-                 porosity_ID
- end enum
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
-   porosity_phasefield_outputID                                                                  !< ID of each post result output
-
-
- public :: &
-   porosity_phasefield_init, &
-   porosity_phasefield_getFormationEnergy, &
-   porosity_phasefield_getSurfaceEnergy, &
-   porosity_phasefield_getSourceAndItsTangent, &
-   porosity_phasefield_getDiffusion33, &
-   porosity_phasefield_getMobility, &
-   porosity_phasefield_putPorosity, &
-   porosity_phasefield_postResults
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine porosity_phasefield_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   porosity_type, &
-   porosity_typeInstance, &
-   homogenization_Noutput, &
-   POROSITY_phasefield_label, &
-   POROSITY_phasefield_ID, &
-   material_homog, & 
-   mappingHomogenization, & 
-   porosityState, &
-   porosityMapping, &
-   porosity, &
-   porosity_initialPhi
- use config, only: &
-   material_partHomogenization, &
-   material_partPhase
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o
- integer(pInt) :: sizeState
- integer(pInt) :: NofMyHomog   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  porosity_'//POROSITY_phasefield_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(porosity_type == POROSITY_phasefield_ID),pInt)
- if (maxNinstance == 0_pInt) return
- 
- allocate(porosity_phasefield_sizePostResults(maxNinstance),                               source=0_pInt)
- allocate(porosity_phasefield_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
- allocate(porosity_phasefield_output         (maxval(homogenization_Noutput),maxNinstance))
-          porosity_phasefield_output = ''
- allocate(porosity_phasefield_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
- allocate(porosity_phasefield_Noutput        (maxNinstance),                               source=0_pInt) 
-
- rewind(fileUnit)
- section = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingHomog: do while (trim(line) /= IO_EOF)                                                       ! read through sections of homog part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next homog section
-     section = section + 1_pInt                                                                     ! advance homog section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (section > 0_pInt ) then; if (porosity_type(section) == POROSITY_phasefield_ID) then          ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = porosity_typeInstance(section)                                                      ! which instance of my porosity is present homog
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('porosity')
-             porosity_phasefield_Noutput(instance) = porosity_phasefield_Noutput(instance) + 1_pInt
-             porosity_phasefield_outputID(porosity_phasefield_Noutput(instance),instance) = porosity_ID
-             porosity_phasefield_output(porosity_phasefield_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
-
-     end select
-   endif; endif
- enddo parsingHomog
- 
- initializeInstances: do section = 1_pInt, size(porosity_type)
-   if (porosity_type(section) == POROSITY_phasefield_ID) then
-     NofMyHomog=count(material_homog==section)
-     instance = porosity_typeInstance(section)
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,porosity_phasefield_Noutput(instance)
-       select case(porosity_phasefield_outputID(o,instance))
-         case(porosity_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          porosity_phasefield_sizePostResult(o,instance) = mySize
-          porosity_phasefield_sizePostResults(instance)  = porosity_phasefield_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-
-! allocate state arrays
-     sizeState = 0_pInt
-     porosityState(section)%sizeState = sizeState
-     porosityState(section)%sizePostResults = porosity_phasefield_sizePostResults(instance)
-     allocate(porosityState(section)%state0   (sizeState,NofMyHomog))
-     allocate(porosityState(section)%subState0(sizeState,NofMyHomog))
-     allocate(porosityState(section)%state    (sizeState,NofMyHomog))
-
-     nullify(porosityMapping(section)%p)
-     porosityMapping(section)%p => mappingHomogenization(1,:,:)
-     deallocate(porosity(section)%p)
-     allocate(porosity(section)%p(NofMyHomog), source=porosity_initialPhi(section))
-     
-   endif
- 
- enddo initializeInstances
-end subroutine porosity_phasefield_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized vacancy formation energy
-!--------------------------------------------------------------------------------------------------
-function porosity_phasefield_getFormationEnergy(ip,el)
- use lattice, only: &
-   lattice_vacancyFormationEnergy, &
-   lattice_vacancyVol
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal) :: &
-   porosity_phasefield_getFormationEnergy
- integer(pInt) :: &
-   grain
-  
- porosity_phasefield_getFormationEnergy = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   porosity_phasefield_getFormationEnergy = porosity_phasefield_getFormationEnergy + &
-    lattice_vacancyFormationEnergy(material_phase(grain,ip,el))/ &
-    lattice_vacancyVol(material_phase(grain,ip,el))
- enddo
-
- porosity_phasefield_getFormationEnergy = &
-   porosity_phasefield_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function porosity_phasefield_getFormationEnergy
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized pore surface energy (normalized by characteristic length)
-!--------------------------------------------------------------------------------------------------
-function porosity_phasefield_getSurfaceEnergy(ip,el)
- use lattice, only: &
-   lattice_vacancySurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal) :: &
-   porosity_phasefield_getSurfaceEnergy
- integer(pInt) :: &
-   grain
-  
- porosity_phasefield_getSurfaceEnergy = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   porosity_phasefield_getSurfaceEnergy = porosity_phasefield_getSurfaceEnergy + &
-    lattice_vacancySurfaceEnergy(material_phase(grain,ip,el))
- enddo
-
- porosity_phasefield_getSurfaceEnergy = &
-   porosity_phasefield_getSurfaceEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function porosity_phasefield_getSurfaceEnergy
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief  calculates homogenized local driving force for pore nucleation and growth  
-!--------------------------------------------------------------------------------------------------
-subroutine porosity_phasefield_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
- use math, only : &
-   math_mul33x33, &
-   math_mul66x6, &
-   math_Mandel33to6, &
-   math_transpose33, &
-   math_I3
- use material, only: &
-   homogenization_Ngrains, &
-   material_homog, &
-   material_phase, &
-   phase_NstiffnessDegradations, &
-   phase_stiffnessDegradation, &
-   vacancyConc, &
-   vacancyfluxMapping, &
-   damage, &
-   damageMapping, &
-   STIFFNESS_DEGRADATION_damage_ID
- use crystallite, only: &
-   crystallite_Fe
- use constitutive, only: &
-   constitutive_homogenizedC  
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   phi
- integer(pInt) :: &
-   phase, &
-   grain, &
-   homog, &
-   mech
- real(pReal) :: &
-   phiDot, dPhiDot_dPhi, Cv, W_e, strain(6), C(6,6)  
-
- homog = material_homog(ip,el)
- Cv = vacancyConc(homog)%p(vacancyfluxMapping(homog)%p(ip,el))
-
- W_e = 0.0_pReal
- do grain = 1, homogenization_Ngrains(homog)
-   phase = material_phase(grain,ip,el)
-   strain = math_Mandel33to6(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,grain,ip,el)), &
-                                           crystallite_Fe(1:3,1:3,grain,ip,el)) - math_I3)/2.0_pReal
-   C = constitutive_homogenizedC(grain,ip,el)
-   do mech = 1_pInt, phase_NstiffnessDegradations(phase)
-     select case(phase_stiffnessDegradation(mech,phase))
-       case (STIFFNESS_DEGRADATION_damage_ID)
-         C = damage(homog)%p(damageMapping(homog)%p(ip,el))* &
-             damage(homog)%p(damageMapping(homog)%p(ip,el))* &
-             C                                        
-     
-     end select
-   enddo                                                 
-   W_e = W_e + sum(abs(strain*math_mul66x6(C,strain)))
- enddo
- W_e = W_e/real(homogenization_Ngrains(homog),pReal)
- 
- phiDot = 2.0_pReal*(1.0_pReal - phi)*(1.0_pReal - Cv)*(1.0_pReal - Cv) - &
-          2.0_pReal*phi*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ &
-                                  porosity_phasefield_getSurfaceEnergy  (ip,el)
- dPhiDot_dPhi = - 2.0_pReal*(1.0_pReal - Cv)*(1.0_pReal - Cv) &
-                - 2.0_pReal*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ &
-                                      porosity_phasefield_getSurfaceEnergy  (ip,el)
-
-end subroutine porosity_phasefield_getSourceAndItsTangent
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized nonlocal diffusion tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function porosity_phasefield_getDiffusion33(ip,el)
- use lattice, only: &
-   lattice_PorosityDiffusion33
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase, &
-   mappingHomogenization
- use crystallite, only: &
-   crystallite_push33ToRef
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   porosity_phasefield_getDiffusion33
- integer(pInt) :: &
-   homog, &
-   grain
-   
- homog  = mappingHomogenization(2,ip,el)
- porosity_phasefield_getDiffusion33 = 0.0_pReal  
- do grain = 1, homogenization_Ngrains(homog)
-   porosity_phasefield_getDiffusion33 = porosity_phasefield_getDiffusion33 + &
-     crystallite_push33ToRef(grain,ip,el,lattice_PorosityDiffusion33(1:3,1:3,material_phase(grain,ip,el)))
- enddo
-
- porosity_phasefield_getDiffusion33 = &
-   porosity_phasefield_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
- 
-end function porosity_phasefield_getDiffusion33
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief Returns homogenized phase field mobility 
-!--------------------------------------------------------------------------------------------------
-real(pReal) function porosity_phasefield_getMobility(ip,el)
- use mesh, only: &
-   mesh_element
- use lattice, only: &
-   lattice_PorosityMobility
- use material, only: &
-   material_phase, &
-   homogenization_Ngrains
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- integer(pInt) :: &
-   ipc
- 
- porosity_phasefield_getMobility = 0.0_pReal
-                                                
- do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
-   porosity_phasefield_getMobility = porosity_phasefield_getMobility &
-                                   + lattice_PorosityMobility(material_phase(ipc,ip,el))
- enddo
-
- porosity_phasefield_getMobility = &
-   porosity_phasefield_getMobility/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
-
-end function porosity_phasefield_getMobility
-
-!--------------------------------------------------------------------------------------------------
-!> @brief updates porosity with solution from phasefield PDE
-!--------------------------------------------------------------------------------------------------
-subroutine porosity_phasefield_putPorosity(phi,ip,el)
- use material, only: &
-   material_homog, &
-   porosityMapping, &
-   porosity
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   phi
- integer(pInt) :: &
-   homog, &
-   offset
- 
- homog  = material_homog(ip,el)
- offset = porosityMapping(homog)%p(ip,el)
- porosity(homog)%p(offset) = phi
-
-end subroutine porosity_phasefield_putPorosity
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of porosity results
-!--------------------------------------------------------------------------------------------------
-function porosity_phasefield_postResults(ip,el)
- use material, only: &
-   mappingHomogenization, &
-   porosity_typeInstance, &
-   porosity
-
- implicit none
- integer(pInt),              intent(in) :: &
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal), dimension(porosity_phasefield_sizePostResults(porosity_typeInstance(mappingHomogenization(2,ip,el)))) :: &
-   porosity_phasefield_postResults
-
- integer(pInt) :: &
-   instance, homog, offset, o, c
-   
- homog     = mappingHomogenization(2,ip,el)
- offset    = mappingHomogenization(1,ip,el)
- instance  = porosity_typeInstance(homog)
-
- c = 0_pInt
- porosity_phasefield_postResults = 0.0_pReal
-
- do o = 1_pInt,porosity_phasefield_Noutput(instance)
-    select case(porosity_phasefield_outputID(o,instance))
- 
-      case (porosity_ID)
-        porosity_phasefield_postResults(c+1_pInt) = porosity(homog)%p(offset)
-        c = c + 1
-    end select
- enddo
-end function porosity_phasefield_postResults
-
-end module porosity_phasefield
--- a/src/source_damage_isoBrittle.f90
+++ b/src/source_damage_isoBrittle.f90
@ -246,10 +246,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
   sourceState, &
   material_homog, &
   phase_NstiffnessDegradations, &
-   phase_stiffnessDegradation, &
-   porosity, &
-   porosityMapping, &  
-   STIFFNESS_DEGRADATION_porosity_ID
+   phase_stiffnessDegradation
 use math, only : &
   math_mul33x33, &
   math_mul66x6, &
@ -279,15 +276,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
 instance = source_damage_isoBrittle_instance(phase)                                                 !< instance of damage_isoBrittle source
 sourceOffset = source_damage_isoBrittle_offset(phase)

- stiffness = C
- do mech = 1_pInt, phase_NstiffnessDegradations(phase)
-   select case(phase_stiffnessDegradation(mech,phase))
-     case (STIFFNESS_DEGRADATION_porosity_ID)
-       stiffness = porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* &
-                   porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* &
-                   stiffness 
-   end select
- enddo                                                 
+ stiffness = C                                            
 strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)

 strainenergy = 2.0_pReal*sum(strain*math_mul66x6(stiffness,strain))/ &
--- a/src/source_vacancy_irradiation.f90
+++ b/src/source_vacancy_irradiation.f90
@ -1,248 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for vacancy generation due to irradiation
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module source_vacancy_irradiation
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_vacancy_irradiation_sizePostResults, &                                                        !< cumulative size of post results
-   source_vacancy_irradiation_offset, &                                                                 !< which source is my current damage mechanism?
-   source_vacancy_irradiation_instance                                                                  !< instance of damage source mechanism
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   source_vacancy_irradiation_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   source_vacancy_irradiation_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   source_vacancy_irradiation_Noutput                                                                   !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),           allocatable,        private :: &
-   source_vacancy_irradiation_cascadeProb, &
-   source_vacancy_irradiation_cascadeVolume
-
- public :: &
-   source_vacancy_irradiation_init, &
-   source_vacancy_irradiation_deltaState, &
-   source_vacancy_irradiation_getRateAndItsTangent
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_irradiation_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use debug, only: &
-   debug_level,&
-   debug_constitutive,&
-   debug_levelBasic
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   phase_source, &
-   phase_Nsources, &
-   phase_Noutput, &
-   SOURCE_vacancy_irradiation_label, &
-   SOURCE_vacancy_irradiation_ID, &
-   material_phase, &  
-   sourceState
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
- use numerics,only: &
-   numerics_integrator
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
- integer(pInt) :: NofMyPhase   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_vacancy_irradiation_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(phase_source == SOURCE_vacancy_irradiation_ID),pInt)
- if (maxNinstance == 0_pInt) return
- if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
- 
- allocate(source_vacancy_irradiation_offset(material_Nphase), source=0_pInt)
- allocate(source_vacancy_irradiation_instance(material_Nphase), source=0_pInt)
- do phase = 1, material_Nphase
-   source_vacancy_irradiation_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_irradiation_ID)
-   do source = 1, phase_Nsources(phase)
-     if (phase_source(source,phase) == source_vacancy_irradiation_ID) &
-       source_vacancy_irradiation_offset(phase) = source
-   enddo    
- enddo
-   
- allocate(source_vacancy_irradiation_sizePostResults(maxNinstance),                     source=0_pInt)
- allocate(source_vacancy_irradiation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
- allocate(source_vacancy_irradiation_output(maxval(phase_Noutput),maxNinstance))
-          source_vacancy_irradiation_output = ''
- allocate(source_vacancy_irradiation_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(source_vacancy_irradiation_cascadeProb(maxNinstance),                         source=0.0_pReal) 
- allocate(source_vacancy_irradiation_cascadeVolume(maxNinstance),                       source=0.0_pReal) 
-
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = source_vacancy_irradiation_instance(phase)                                          ! which instance of my vacancy is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('irradiation_cascadeprobability')
-         source_vacancy_irradiation_cascadeProb(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-       case ('irradiation_cascadevolume')
-         source_vacancy_irradiation_cascadeVolume(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-     end select
-   endif; endif
- enddo parsingFile
- 
- initializeInstances: do phase = 1_pInt, material_Nphase
-   if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then
-     NofMyPhase=count(material_phase==phase)
-     instance = source_vacancy_irradiation_instance(phase)
-     sourceOffset = source_vacancy_irradiation_offset(phase)
-
-     sizeDotState              =   2_pInt
-     sizeDeltaState            =   2_pInt
-     sizeState                 =   2_pInt
-     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_irradiation_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.1_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
-
-   endif
- 
- enddo initializeInstances
-end subroutine source_vacancy_irradiation_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_irradiation_deltaState(ipc, ip, el)
- use material, only: &
-   phaseAt, phasememberAt, &
-   sourceState
-
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- integer(pInt) :: &
-   phase, constituent, sourceOffset
- real(pReal) :: &
-   randNo   
-
- phase       = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
- sourceOffset = source_vacancy_irradiation_offset(phase)
-
- call random_number(randNo)
- sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
-   randNo - sourceState(phase)%p(sourceOffset)%state(1,constituent)
- call random_number(randNo)
- sourceState(phase)%p(sourceOffset)%deltaState(2,constituent) = &
-   randNo - sourceState(phase)%p(sourceOffset)%state(2,constituent)
-
-end subroutine source_vacancy_irradiation_deltaState
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local vacancy generation rate 
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_irradiation_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
- use material, only: &
-   phaseAt, phasememberAt, &
-   sourceState
-
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), intent(out) :: &
-   CvDot, dCvDot_dCv
- integer(pInt) :: &
-   instance, phase, constituent, sourceOffset 
-
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
- instance = source_vacancy_irradiation_instance(phase)
- sourceOffset = source_vacancy_irradiation_offset(phase)
- 
- CvDot = 0.0_pReal
- dCvDot_dCv = 0.0_pReal
- if (sourceState(phase)%p(sourceOffset)%state0(1,constituent) < source_vacancy_irradiation_cascadeProb(instance)) &
-   CvDot = sourceState(phase)%p(sourceOffset)%state0(2,constituent)*source_vacancy_irradiation_cascadeVolume(instance)
- 
-end subroutine source_vacancy_irradiation_getRateAndItsTangent
- 
-end module source_vacancy_irradiation
--- a/src/source_vacancy_phenoplasticity.f90
+++ b/src/source_vacancy_phenoplasticity.f90
@ -1,210 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for vacancy generation due to plasticity
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module source_vacancy_phenoplasticity
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_vacancy_phenoplasticity_sizePostResults, &                                                        !< cumulative size of post results
-   source_vacancy_phenoplasticity_offset, &                                                                 !< which source is my current damage mechanism?
-   source_vacancy_phenoplasticity_instance                                                                  !< instance of damage source mechanism
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   source_vacancy_phenoplasticity_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   source_vacancy_phenoplasticity_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   source_vacancy_phenoplasticity_Noutput                                                                   !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),           allocatable,        private :: &
-   source_vacancy_phenoplasticity_rateCoeff
-
- public :: &
-   source_vacancy_phenoplasticity_init, &
-   source_vacancy_phenoplasticity_getRateAndItsTangent
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_phenoplasticity_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use debug, only: &
-   debug_level,&
-   debug_constitutive,&
-   debug_levelBasic
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   phase_source, &
-   phase_Nsources, &
-   phase_Noutput, &
-   SOURCE_vacancy_phenoplasticity_label, &
-   SOURCE_vacancy_phenoplasticity_ID, &
-   material_phase, &  
-   sourceState
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
- use numerics,only: &
-   numerics_integrator
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
- integer(pInt) :: NofMyPhase   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_vacancy_phenoplasticity_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(phase_source == SOURCE_vacancy_phenoplasticity_ID),pInt)
- if (maxNinstance == 0_pInt) return
- if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
- 
- allocate(source_vacancy_phenoplasticity_offset(material_Nphase), source=0_pInt)
- allocate(source_vacancy_phenoplasticity_instance(material_Nphase), source=0_pInt)
- do phase = 1, material_Nphase
-   source_vacancy_phenoplasticity_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_phenoplasticity_ID)
-   do source = 1, phase_Nsources(phase)
-     if (phase_source(source,phase) == source_vacancy_phenoplasticity_ID) &
-       source_vacancy_phenoplasticity_offset(phase) = source
-   enddo    
- enddo
-   
- allocate(source_vacancy_phenoplasticity_sizePostResults(maxNinstance),                     source=0_pInt)
- allocate(source_vacancy_phenoplasticity_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
- allocate(source_vacancy_phenoplasticity_output(maxval(phase_Noutput),maxNinstance))
-          source_vacancy_phenoplasticity_output = ''
- allocate(source_vacancy_phenoplasticity_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(source_vacancy_phenoplasticity_rateCoeff(maxNinstance),                           source=0.0_pReal) 
-
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = source_vacancy_phenoplasticity_instance(phase)                                          ! which instance of my vacancy is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('phenoplasticity_ratecoeff')
-         source_vacancy_phenoplasticity_rateCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-     end select
-   endif; endif
- enddo parsingFile
- 
- initializeInstances: do phase = 1_pInt, material_Nphase
-   if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then
-     NofMyPhase=count(material_phase==phase)
-     instance = source_vacancy_phenoplasticity_instance(phase)
-     sourceOffset = source_vacancy_phenoplasticity_offset(phase)
-
-     sizeDotState              =   0_pInt
-     sizeDeltaState            =   0_pInt
-     sizeState                 =   0_pInt
-     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_phenoplasticity_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
-
-   endif
- 
- enddo initializeInstances
-end subroutine source_vacancy_phenoplasticity_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local vacancy generation rate 
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_phenoplasticity_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
- use material, only: &
-   phaseAt, phasememberAt, &
-   plasticState
-
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), intent(out) :: &
-   CvDot, dCvDot_dCv
- integer(pInt) :: &
-   instance, phase, constituent 
-
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
- instance = source_vacancy_phenoplasticity_instance(phase)
- 
- CvDot = &
-   source_vacancy_phenoplasticity_rateCoeff(instance)* &
-   sum(plasticState(phase)%slipRate(:,constituent))
- dCvDot_dCv = 0.0_pReal  
- 
-end subroutine source_vacancy_phenoplasticity_getRateAndItsTangent
- 
-end module source_vacancy_phenoplasticity
--- a/src/source_vacancy_thermalfluc.f90
+++ b/src/source_vacancy_thermalfluc.f90
@ -1,250 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for vacancy generation due to thermal fluctuations
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module source_vacancy_thermalfluc
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_vacancy_thermalfluc_sizePostResults, &                                                        !< cumulative size of post results
-   source_vacancy_thermalfluc_offset, &                                                                 !< which source is my current damage mechanism?
-   source_vacancy_thermalfluc_instance                                                                  !< instance of damage source mechanism
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   source_vacancy_thermalfluc_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   source_vacancy_thermalfluc_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   source_vacancy_thermalfluc_Noutput                                                                   !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),           allocatable,        private :: &
-   source_vacancy_thermalfluc_amplitude, &
-   source_vacancy_thermalfluc_normVacancyEnergy
-
- public :: &
-   source_vacancy_thermalfluc_init, &
-   source_vacancy_thermalfluc_deltaState, &
-   source_vacancy_thermalfluc_getRateAndItsTangent
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_thermalfluc_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use debug, only: &
-   debug_level,&
-   debug_constitutive,&
-   debug_levelBasic
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use lattice, only: &
-   lattice_vacancyFormationEnergy
- use material, only: &
-   phase_source, &
-   phase_Nsources, &
-   phase_Noutput, &
-   SOURCE_vacancy_thermalfluc_label, &
-   SOURCE_vacancy_thermalfluc_ID, &
-   material_phase, &  
-   sourceState
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
- use numerics,only: &
-   numerics_integrator
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
- integer(pInt) :: NofMyPhase   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_vacancy_thermalfluc_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(phase_source == SOURCE_vacancy_thermalfluc_ID),pInt)
- if (maxNinstance == 0_pInt) return
- if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
- 
- allocate(source_vacancy_thermalfluc_offset(material_Nphase), source=0_pInt)
- allocate(source_vacancy_thermalfluc_instance(material_Nphase), source=0_pInt)
- do phase = 1, material_Nphase
-   source_vacancy_thermalfluc_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_thermalfluc_ID)
-   do source = 1, phase_Nsources(phase)
-     if (phase_source(source,phase) == source_vacancy_thermalfluc_ID) &
-       source_vacancy_thermalfluc_offset(phase) = source
-   enddo    
- enddo
-   
- allocate(source_vacancy_thermalfluc_sizePostResults(maxNinstance),                     source=0_pInt)
- allocate(source_vacancy_thermalfluc_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
- allocate(source_vacancy_thermalfluc_output(maxval(phase_Noutput),maxNinstance))
-          source_vacancy_thermalfluc_output = ''
- allocate(source_vacancy_thermalfluc_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(source_vacancy_thermalfluc_amplitude(maxNinstance),                           source=0.0_pReal) 
- allocate(source_vacancy_thermalfluc_normVacancyEnergy(maxNinstance),                   source=0.0_pReal) 
-
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = source_vacancy_thermalfluc_instance(phase)                                          ! which instance of my vacancy is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('thermalfluctuation_amplitude')
-         source_vacancy_thermalfluc_amplitude(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-     end select
-   endif; endif
- enddo parsingFile
- 
- initializeInstances: do phase = 1_pInt, material_Nphase
-   if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then
-     NofMyPhase=count(material_phase==phase)
-     instance = source_vacancy_thermalfluc_instance(phase)
-     source_vacancy_thermalfluc_normVacancyEnergy(instance) = &
-       lattice_vacancyFormationEnergy(phase)/1.3806488e-23_pReal
-     sourceOffset = source_vacancy_thermalfluc_offset(phase)
-
-     sizeDotState              =   1_pInt
-     sizeDeltaState            =   1_pInt
-     sizeState                 =   1_pInt
-     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.1_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
-
-   endif
- 
- enddo initializeInstances
-end subroutine source_vacancy_thermalfluc_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_thermalfluc_deltaState(ipc, ip, el)
- use material, only: &
-   phaseAt, phasememberAt, &
-   sourceState
-
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- integer(pInt) :: &
-   phase, constituent, sourceOffset 
- real(pReal) :: &
-   randNo   
-
- phase       = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
- sourceOffset = source_vacancy_thermalfluc_offset(phase)
-
- call random_number(randNo)
- sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
-   randNo - 0.5_pReal - sourceState(phase)%p(sourceOffset)%state(1,constituent)
-
-end subroutine source_vacancy_thermalfluc_deltaState
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local vacancy generation rate 
-!--------------------------------------------------------------------------------------------------
-subroutine source_vacancy_thermalfluc_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
- use material, only: &
-   phaseAt, phasememberAt, &
-   material_homog, &
-   temperature, &
-   thermalMapping, &
-   sourceState
-
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< grain number
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), intent(out) :: &
-   CvDot, dCvDot_dCv
- integer(pInt) :: &
-   instance, phase, constituent, sourceOffset 
-
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
- instance = source_vacancy_thermalfluc_instance(phase)
- sourceOffset = source_vacancy_thermalfluc_offset(phase)
- 
- CvDot = source_vacancy_thermalfluc_amplitude(instance)* &
-         sourceState(phase)%p(sourceOffset)%state0(2,constituent)* &
-         exp(-source_vacancy_thermalfluc_normVacancyEnergy(instance)/ &
-              temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el)))
- dCvDot_dCv = 0.0_pReal
- 
-end subroutine source_vacancy_thermalfluc_getRateAndItsTangent
- 
-end module source_vacancy_thermalfluc
--- a/src/vacancyflux_cahnhilliard.f90
+++ b/src/vacancyflux_cahnhilliard.f90
@ -1,602 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for conservative transport of vacancy concentration field
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module vacancyflux_cahnhilliard
- use prec, only: &
-   pReal, &
-   pInt, &
-   group_float
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   vacancyflux_cahnhilliard_sizePostResults                                                           !< cumulative size of post results
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   vacancyflux_cahnhilliard_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   vacancyflux_cahnhilliard_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   vacancyflux_cahnhilliard_Noutput                                                                   !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),           allocatable,        private :: &
-   vacancyflux_cahnhilliard_flucAmplitude
-
- type(group_float),                  dimension(:),           allocatable,        private :: &
-   vacancyflux_cahnhilliard_thermalFluc
-
- real(pReal),                                                 parameter,           private :: &
-   kB = 1.3806488e-23_pReal                                                                          !< Boltzmann constant in J/Kelvin
-
- enum, bind(c) 
-   enumerator :: undefined_ID, &
-                 vacancyConc_ID
- end enum
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,        private :: & 
-   vacancyflux_cahnhilliard_outputID                                                                  !< ID of each post result output
-
-
- public :: &
-   vacancyflux_cahnhilliard_init, &
-   vacancyflux_cahnhilliard_getSourceAndItsTangent, &
-   vacancyflux_cahnhilliard_getMobility33, &
-   vacancyflux_cahnhilliard_getDiffusion33, &
-   vacancyflux_cahnhilliard_getChemPotAndItsTangent, &
-   vacancyflux_cahnhilliard_putVacancyConcAndItsRate, &
-   vacancyflux_cahnhilliard_postResults
- private :: &
-   vacancyflux_cahnhilliard_getFormationEnergy, &
-   vacancyflux_cahnhilliard_getEntropicCoeff, &
-   vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent  
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_cahnhilliard_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   vacancyflux_type, &
-   vacancyflux_typeInstance, &
-   homogenization_Noutput, &
-   VACANCYFLUX_cahnhilliard_label, &
-   VACANCYFLUX_cahnhilliard_ID, &
-   material_homog, &  
-   mappingHomogenization, &
-   vacancyfluxState, &
-   vacancyfluxMapping, &
-   vacancyConc, &
-   vacancyConcRate, &
-   vacancyflux_initialCv
- use config, only: &
-   material_partPhase, &
-   material_partHomogenization
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o,offset
- integer(pInt) :: sizeState
- integer(pInt) :: NofMyHomog   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID),pInt)
- if (maxNinstance == 0_pInt) return
- 
- allocate(vacancyflux_cahnhilliard_sizePostResults(maxNinstance),                               source=0_pInt)
- allocate(vacancyflux_cahnhilliard_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
- allocate(vacancyflux_cahnhilliard_output         (maxval(homogenization_Noutput),maxNinstance))
-          vacancyflux_cahnhilliard_output = ''
- allocate(vacancyflux_cahnhilliard_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
- allocate(vacancyflux_cahnhilliard_Noutput        (maxNinstance),                               source=0_pInt) 
-
- allocate(vacancyflux_cahnhilliard_flucAmplitude  (maxNinstance))
- allocate(vacancyflux_cahnhilliard_thermalFluc    (maxNinstance))
-
- rewind(fileUnit)
- section = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingHomog: do while (trim(line) /= IO_EOF)                                                      ! read through sections of homog part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next homog section
-     section = section + 1_pInt                                                                     ! advance homog section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (section > 0_pInt ) then; if (vacancyflux_type(section) == VACANCYFLUX_cahnhilliard_ID) then  ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = vacancyflux_typeInstance(section)                                                   ! which instance of my vacancyflux is present homog
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('vacancyconc')
-             vacancyflux_cahnhilliard_Noutput(instance) = vacancyflux_cahnhilliard_Noutput(instance) + 1_pInt
-             vacancyflux_cahnhilliard_outputID(vacancyflux_cahnhilliard_Noutput(instance),instance) = vacancyConc_ID
-             vacancyflux_cahnhilliard_output(vacancyflux_cahnhilliard_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
-
-       case ('vacancyflux_flucamplitude')
-         vacancyflux_cahnhilliard_flucAmplitude(instance) = IO_floatValue(line,chunkPos,2_pInt)
-         
-     end select
-   endif; endif
- enddo parsingHomog
- 
- initializeInstances: do section = 1_pInt, size(vacancyflux_type)
-   if (vacancyflux_type(section) == VACANCYFLUX_cahnhilliard_ID) then
-     NofMyHomog=count(material_homog==section)
-     instance = vacancyflux_typeInstance(section)
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,vacancyflux_cahnhilliard_Noutput(instance)
-       select case(vacancyflux_cahnhilliard_outputID(o,instance))
-         case(vacancyConc_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          vacancyflux_cahnhilliard_sizePostResult(o,instance) = mySize
-          vacancyflux_cahnhilliard_sizePostResults(instance)  = vacancyflux_cahnhilliard_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-
-! allocate state arrays
-     sizeState = 0_pInt
-     vacancyfluxState(section)%sizeState = sizeState
-     vacancyfluxState(section)%sizePostResults = vacancyflux_cahnhilliard_sizePostResults(instance)
-     allocate(vacancyfluxState(section)%state0   (sizeState,NofMyHomog))
-     allocate(vacancyfluxState(section)%subState0(sizeState,NofMyHomog))
-     allocate(vacancyfluxState(section)%state    (sizeState,NofMyHomog))
-     
-     allocate(vacancyflux_cahnhilliard_thermalFluc(instance)%p(NofMyHomog))
-     do offset = 1_pInt, NofMyHomog
-       call random_number(vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset))
-       vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset) = &
-         1.0_pReal - &
-         vacancyflux_cahnhilliard_flucAmplitude(instance)* &
-         (vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset) - 0.5_pReal)
-     enddo  
-
-     nullify(vacancyfluxMapping(section)%p)
-     vacancyfluxMapping(section)%p => mappingHomogenization(1,:,:)
-     deallocate(vacancyConc    (section)%p)
-     allocate  (vacancyConc    (section)%p(NofMyHomog), source=vacancyflux_initialCv(section))
-     deallocate(vacancyConcRate(section)%p)
-     allocate  (vacancyConcRate(section)%p(NofMyHomog), source=0.0_pReal)
-     
-   endif
- 
- enddo initializeInstances
- 
-end subroutine vacancyflux_cahnhilliard_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  calculates homogenized vacancy driving forces  
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el)
- use material, only: &
-   homogenization_Ngrains, &
-   mappingHomogenization, &
-   phaseAt, &
-   phase_source, &
-   phase_Nsources, &
-   SOURCE_vacancy_phenoplasticity_ID, &
-   SOURCE_vacancy_irradiation_ID, &
-   SOURCE_vacancy_thermalfluc_ID
- use source_vacancy_phenoplasticity, only: &
-   source_vacancy_phenoplasticity_getRateAndItsTangent
- use source_vacancy_irradiation, only: &
-   source_vacancy_irradiation_getRateAndItsTangent
- use source_vacancy_thermalfluc, only: &
-   source_vacancy_thermalfluc_getRateAndItsTangent
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   Cv
- integer(pInt) :: &
-   phase, &
-   grain, &
-   source
- real(pReal) :: &
-   CvDot, dCvDot_dCv, localCvDot, dLocalCvDot_dCv  
-
- CvDot = 0.0_pReal
- dCvDot_dCv = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el))
-   phase = phaseAt(grain,ip,el)
-   do source = 1_pInt, phase_Nsources(phase)
-     select case(phase_source(source,phase))                                                   
-       case (SOURCE_vacancy_phenoplasticity_ID)
-        call source_vacancy_phenoplasticity_getRateAndItsTangent  (localCvDot, dLocalCvDot_dCv, grain, ip,  el)
-
-       case (SOURCE_vacancy_irradiation_ID)
-        call source_vacancy_irradiation_getRateAndItsTangent  (localCvDot, dLocalCvDot_dCv, grain, ip,  el)
-
-       case (SOURCE_vacancy_thermalfluc_ID)
-        call source_vacancy_thermalfluc_getRateAndItsTangent(localCvDot, dLocalCvDot_dCv, grain, ip,  el)
-
-     end select
-     CvDot = CvDot + localCvDot
-     dCvDot_dCv = dCvDot_dCv + dLocalCvDot_dCv
-   enddo  
- enddo
- 
- CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
- dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
- 
-end subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized vacancy mobility tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function vacancyflux_cahnhilliard_getMobility33(ip,el)
- use lattice, only: &
-   lattice_vacancyfluxMobility33
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_push33ToRef
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   vacancyflux_cahnhilliard_getMobility33
- integer(pInt) :: &
-   grain
-  
- vacancyflux_cahnhilliard_getMobility33 = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   vacancyflux_cahnhilliard_getMobility33 = vacancyflux_cahnhilliard_getMobility33 + &
-    crystallite_push33ToRef(grain,ip,el,lattice_vacancyfluxMobility33(:,:,material_phase(grain,ip,el)))
- enddo
-
- vacancyflux_cahnhilliard_getMobility33 = &
-   vacancyflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function vacancyflux_cahnhilliard_getMobility33
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized vacancy diffusion tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function vacancyflux_cahnhilliard_getDiffusion33(ip,el)
- use lattice, only: &
-   lattice_vacancyfluxDiffusion33
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_push33ToRef
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal), dimension(3,3) :: &
-   vacancyflux_cahnhilliard_getDiffusion33
- integer(pInt) :: &
-   grain
-  
- vacancyflux_cahnhilliard_getDiffusion33 = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   vacancyflux_cahnhilliard_getDiffusion33 = vacancyflux_cahnhilliard_getDiffusion33 + &
-    crystallite_push33ToRef(grain,ip,el,lattice_vacancyfluxDiffusion33(:,:,material_phase(grain,ip,el)))
- enddo
-
- vacancyflux_cahnhilliard_getDiffusion33 = &
-   vacancyflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function vacancyflux_cahnhilliard_getDiffusion33
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized vacancy formation energy
-!--------------------------------------------------------------------------------------------------
-real(pReal) function vacancyflux_cahnhilliard_getFormationEnergy(ip,el)
- use lattice, only: &
-   lattice_vacancyFormationEnergy, &
-   lattice_vacancyVol, &
-   lattice_vacancySurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_phase
- use mesh, only: &
-   mesh_element
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- integer(pInt) :: &
-   grain
-  
- vacancyflux_cahnhilliard_getFormationEnergy = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mesh_element(3,el))
-   vacancyflux_cahnhilliard_getFormationEnergy = vacancyflux_cahnhilliard_getFormationEnergy + &
-    lattice_vacancyFormationEnergy(material_phase(grain,ip,el))/ &
-    lattice_vacancyVol(material_phase(grain,ip,el))/ &
-    lattice_vacancySurfaceEnergy(material_phase(grain,ip,el))
- enddo
-
- vacancyflux_cahnhilliard_getFormationEnergy = &
-   vacancyflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
- 
-end function vacancyflux_cahnhilliard_getFormationEnergy
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized vacancy entropy coefficient
-!--------------------------------------------------------------------------------------------------
-real(pReal) function vacancyflux_cahnhilliard_getEntropicCoeff(ip,el)
- use lattice, only: &
-   lattice_vacancyVol, &
-   lattice_vacancySurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_homog, &
-   material_phase, &
-   temperature, &
-   thermalMapping
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- integer(pInt) :: &
-   grain
-  
- vacancyflux_cahnhilliard_getEntropicCoeff = 0.0_pReal
- do grain = 1, homogenization_Ngrains(material_homog(ip,el))
-   vacancyflux_cahnhilliard_getEntropicCoeff = vacancyflux_cahnhilliard_getEntropicCoeff + &
-    kB/ &
-    lattice_vacancyVol(material_phase(grain,ip,el))/ &
-    lattice_vacancySurfaceEnergy(material_phase(grain,ip,el))
- enddo
-
- vacancyflux_cahnhilliard_getEntropicCoeff = &
-   vacancyflux_cahnhilliard_getEntropicCoeff* &
-   temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ &
-   real(homogenization_Ngrains(material_homog(ip,el)),pReal)
- 
-end function vacancyflux_cahnhilliard_getEntropicCoeff
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized kinematic contribution to chemical potential
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCv, Cv, ip, el)
- use lattice, only: &
-   lattice_vacancySurfaceEnergy
- use material, only: &
-   homogenization_Ngrains, &
-   material_homog, &
-   phase_kinematics, &
-   phase_Nkinematics, &
-   material_phase, &
-   KINEMATICS_vacancy_strain_ID
- use crystallite, only: &
-   crystallite_Tstar_v, &
-   crystallite_Fi0, &
-   crystallite_Fi
- use kinematics_vacancy_strain, only: &
-   kinematics_vacancy_strain_ChemPotAndItsTangent
-
- implicit none
- integer(pInt), intent(in)  :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in)  :: &
-   Cv
- real(pReal),   intent(out) :: &
-   KPot, dKPot_dCv
- real(pReal) :: &
-   my_KPot, my_dKPot_dCv
- integer(pInt) :: &
-   grain, kinematics
-  
- KPot = 0.0_pReal
- dKPot_dCv = 0.0_pReal
- do grain = 1_pInt,homogenization_Ngrains(material_homog(ip,el)) 
-   do kinematics = 1_pInt, phase_Nkinematics(material_phase(grain,ip,el))
-     select case (phase_kinematics(kinematics,material_phase(grain,ip,el)))
-       case (KINEMATICS_vacancy_strain_ID)
-         call kinematics_vacancy_strain_ChemPotAndItsTangent(my_KPot, my_dKPot_dCv, &
-                                                             crystallite_Tstar_v(1:6,grain,ip,el), &
-                                                             crystallite_Fi0(1:3,1:3,grain,ip,el), &
-                                                             crystallite_Fi (1:3,1:3,grain,ip,el), &
-                                                             grain,ip, el)
-
-       case default
-         my_KPot = 0.0_pReal
-         my_dKPot_dCv = 0.0_pReal
-     
-     end select
-     KPot = KPot + my_KPot/lattice_vacancySurfaceEnergy(material_phase(grain,ip,el))
-     dKPot_dCv = dKPot_dCv + my_dKPot_dCv/lattice_vacancySurfaceEnergy(material_phase(grain,ip,el))
-   enddo
- enddo 
- 
- KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
- dKPot_dCv = dKPot_dCv/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
-
-end subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized chemical potential and its tangent
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_cahnhilliard_getChemPotAndItsTangent(ChemPot,dChemPot_dCv,Cv,ip,el)
- use numerics, only: &
-   vacancyBoundPenalty, &
-   vacancyPolyOrder
- use material, only: &
-   mappingHomogenization, &
-   vacancyflux_typeInstance, &
-   porosity, &
-   porosityMapping
-
- implicit none
- integer(pInt), intent(in)  :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in)  :: &
-   Cv
- real(pReal),   intent(out) :: &
-   ChemPot, &
-   dChemPot_dCv
- real(pReal) :: &
-   VoidPhaseFrac, kBT, KPot, dKPot_dCv  
- integer(pInt) :: &
-   homog, o  
-  
- homog = mappingHomogenization(2,ip,el)
- VoidPhaseFrac = porosity(homog)%p(porosityMapping(homog)%p(ip,el))
- kBT = vacancyflux_cahnhilliard_getEntropicCoeff(ip,el)
- 
- ChemPot = vacancyflux_cahnhilliard_getFormationEnergy(ip,el)    
- dChemPot_dCv = 0.0_pReal
- do o = 1_pInt, vacancyPolyOrder
-   ChemPot = ChemPot + kBT*((2.0_pReal*Cv - 1.0_pReal)**real(2_pInt*o-1_pInt,pReal))/ &
-                       real(2_pInt*o-1_pInt,pReal)
-   dChemPot_dCv = dChemPot_dCv + 2.0_pReal*kBT*(2.0_pReal*Cv - 1.0_pReal)**real(2_pInt*o-2_pInt,pReal) 
- enddo
- 
- ChemPot =   VoidPhaseFrac*VoidPhaseFrac*ChemPot &
-           - 2.0_pReal*(1.0_pReal - Cv)*(1.0_pReal - VoidPhaseFrac)*(1.0_pReal - VoidPhaseFrac)
-
- dChemPot_dCv =   VoidPhaseFrac*VoidPhaseFrac*dChemPot_dCv &
-                + 2.0_pReal*(1.0_pReal - VoidPhaseFrac)*(1.0_pReal - VoidPhaseFrac)
-
- call vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCv, Cv, ip, el)
- ChemPot = ChemPot + KPot 
- dChemPot_dCv = dChemPot_dCv + dKPot_dCv
- 
- if (Cv < 0.0_pReal) then
-   ChemPot = ChemPot - 3.0_pReal*vacancyBoundPenalty*Cv*Cv
-   dChemPot_dCv = dChemPot_dCv - 6.0_pReal*vacancyBoundPenalty*Cv
- elseif (Cv > 1.0_pReal) then
-   ChemPot = ChemPot + 3.0_pReal*vacancyBoundPenalty*(1.0_pReal - Cv)*(1.0_pReal - Cv)
-   dChemPot_dCv = dChemPot_dCv - 6.0_pReal*vacancyBoundPenalty*(1.0_pReal - Cv)
- endif
- 
- ChemPot = ChemPot* &
-           vacancyflux_cahnhilliard_thermalFluc(vacancyflux_typeInstance(homog))%p(mappingHomogenization(1,ip,el))
- dChemPot_dCv = dChemPot_dCv* &
-           vacancyflux_cahnhilliard_thermalFluc(vacancyflux_typeInstance(homog))%p(mappingHomogenization(1,ip,el))                 
- 
-end subroutine vacancyflux_cahnhilliard_getChemPotAndItsTangent
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief updated vacancy concentration and its rate with solution from transport PDE
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_cahnhilliard_putVacancyConcAndItsRate(Cv,Cvdot,ip,el)
- use material, only: &
-   mappingHomogenization, &
-   vacancyConc, &
-   vacancyConcRate, &
-   vacancyfluxMapping
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   Cv, &
-   Cvdot
- integer(pInt) :: &
-   homog, &
-   offset  
- 
- homog  = mappingHomogenization(2,ip,el)
- offset = vacancyfluxMapping(homog)%p(ip,el)
- vacancyConc    (homog)%p(offset) = Cv
- vacancyConcRate(homog)%p(offset) = Cvdot
-
-end subroutine vacancyflux_cahnhilliard_putVacancyConcAndItsRate
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of vacancy transport results
-!--------------------------------------------------------------------------------------------------
-function vacancyflux_cahnhilliard_postResults(ip,el)
- use material, only: &
-   mappingHomogenization, &
-   vacancyflux_typeInstance, &
-   vacancyConc, &
-   vacancyfluxMapping
-
- implicit none
- integer(pInt),              intent(in) :: &
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal), dimension(vacancyflux_cahnhilliard_sizePostResults(vacancyflux_typeInstance(mappingHomogenization(2,ip,el)))) :: &
-   vacancyflux_cahnhilliard_postResults
-
- integer(pInt) :: &
-   instance, homog, offset, o, c
-   
- homog     = mappingHomogenization(2,ip,el)
- offset    = vacancyfluxMapping(homog)%p(ip,el)
- instance  = vacancyflux_typeInstance(homog)
-
- c = 0_pInt
- vacancyflux_cahnhilliard_postResults = 0.0_pReal
-
- do o = 1_pInt,vacancyflux_cahnhilliard_Noutput(instance)
-    select case(vacancyflux_cahnhilliard_outputID(o,instance))
- 
-      case (vacancyConc_ID)
-        vacancyflux_cahnhilliard_postResults(c+1_pInt) = vacancyConc(homog)%p(offset)
-        c = c + 1
-    end select
- enddo
-end function vacancyflux_cahnhilliard_postResults
-
-end module vacancyflux_cahnhilliard
--- a/src/vacancyflux_isochempot.f90
+++ b/src/vacancyflux_isochempot.f90
@ -1,328 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for locally evolving vacancy concentration
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-module vacancyflux_isochempot
- use prec, only: &
-   pReal, &
-   pInt
-
- implicit none
- private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   vacancyflux_isochempot_sizePostResults                                                           !< cumulative size of post results
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
-   vacancyflux_isochempot_sizePostResult                                                            !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public :: &
-   vacancyflux_isochempot_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   vacancyflux_isochempot_Noutput                                                                   !< number of outputs per instance of this damage 
-
- enum, bind(c) 
-   enumerator :: undefined_ID, &
-                 vacancyconc_ID
- end enum
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
-   vacancyflux_isochempot_outputID                                                                  !< ID of each post result output
-
-
- public :: &
-   vacancyflux_isochempot_init, &
-   vacancyflux_isochempot_updateState, &
-   vacancyflux_isochempot_getSourceAndItsTangent, &
-   vacancyflux_isochempot_postResults 
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_isochempot_init(fileUnit)
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
- use material, only: &
-   vacancyflux_type, &
-   vacancyflux_typeInstance, &
-   homogenization_Noutput, &
-   VACANCYFLUX_isochempot_label, &
-   VACANCYFLUX_isochempot_ID, &
-   material_homog, &  
-   mappingHomogenization, &
-   vacancyfluxState, &
-   vacancyfluxMapping, &
-   vacancyConc, &
-   vacancyConcRate, &
-   vacancyflux_initialCv
- use config, only: &
-   material_partHomogenization
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o
- integer(pInt) :: sizeState
- integer(pInt) :: NofMyHomog   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-
- write(6,'(/,a)')   ' <<<+-  vacancyflux_'//VACANCYFLUX_isochempot_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- 
- maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_isochempot_ID),pInt)
- if (maxNinstance == 0_pInt) return
- 
- allocate(vacancyflux_isochempot_sizePostResults(maxNinstance),                               source=0_pInt)
- allocate(vacancyflux_isochempot_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
- allocate(vacancyflux_isochempot_output         (maxval(homogenization_Noutput),maxNinstance))
-          vacancyflux_isochempot_output = ''
- allocate(vacancyflux_isochempot_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
- allocate(vacancyflux_isochempot_Noutput        (maxNinstance),                               source=0_pInt) 
-
- rewind(fileUnit)
- section = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of homog part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next homog section
-     section = section + 1_pInt                                                                     ! advance homog section counter
-     cycle                                                                                          ! skip to next line
-   endif
-
-   if (section > 0_pInt ) then; if (vacancyflux_type(section) == VACANCYFLUX_isochempot_ID) then             ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-
-     instance = vacancyflux_typeInstance(section)                                                       ! which instance of my vacancyflux is present homog
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('vacancyconc')
-             vacancyflux_isochempot_Noutput(instance) = vacancyflux_isochempot_Noutput(instance) + 1_pInt
-             vacancyflux_isochempot_outputID(vacancyflux_isochempot_Noutput(instance),instance) = vacancyconc_ID
-             vacancyflux_isochempot_output(vacancyflux_isochempot_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
-
-     end select
-   endif; endif
- enddo parsingFile
- 
- initializeInstances: do section = 1_pInt, size(vacancyflux_type)
-   if (vacancyflux_type(section) == VACANCYFLUX_isochempot_ID) then
-     NofMyHomog=count(material_homog==section)
-     instance = vacancyflux_typeInstance(section)
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,vacancyflux_isochempot_Noutput(instance)
-       select case(vacancyflux_isochempot_outputID(o,instance))
-         case(vacancyconc_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          vacancyflux_isochempot_sizePostResult(o,instance) = mySize
-          vacancyflux_isochempot_sizePostResults(instance)  = vacancyflux_isochempot_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-
-! allocate state arrays
-     sizeState = 1_pInt
-     vacancyfluxState(section)%sizeState = sizeState
-     vacancyfluxState(section)%sizePostResults = vacancyflux_isochempot_sizePostResults(instance)
-     allocate(vacancyfluxState(section)%state0   (sizeState,NofMyHomog), source=vacancyflux_initialCv(section))
-     allocate(vacancyfluxState(section)%subState0(sizeState,NofMyHomog), source=vacancyflux_initialCv(section))
-     allocate(vacancyfluxState(section)%state    (sizeState,NofMyHomog), source=vacancyflux_initialCv(section))
-
-     nullify(vacancyfluxMapping(section)%p)
-     vacancyfluxMapping(section)%p => mappingHomogenization(1,:,:)
-     deallocate(vacancyConc(section)%p)
-     vacancyConc(section)%p => vacancyfluxState(section)%state(1,:)
-     deallocate(vacancyConcRate(section)%p)
-     allocate(vacancyConcRate(section)%p(NofMyHomog), source=0.0_pReal)
-     
-   endif
- 
- enddo initializeInstances
-end subroutine vacancyflux_isochempot_init
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  calculates change in vacancy concentration based on local vacancy generation model  
-!--------------------------------------------------------------------------------------------------
-function vacancyflux_isochempot_updateState(subdt, ip, el)
- use numerics, only: &
-   err_vacancyflux_tolAbs, &
-   err_vacancyflux_tolRel
- use material, only: &
-   mappingHomogenization, &
-   vacancyflux_typeInstance, &
-   vacancyfluxState, &
-   vacancyConc, &
-   vacancyConcRate, &
-   vacancyfluxMapping
-
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   subdt
- logical,                    dimension(2)                             :: &
-   vacancyflux_isochempot_updateState
- integer(pInt) :: &
-   homog, &
-   offset, &
-   instance
- real(pReal) :: &
-   Cv, Cvdot, dCvDot_dCv  
-
- homog  = mappingHomogenization(2,ip,el)
- offset = mappingHomogenization(1,ip,el)
- instance = vacancyflux_typeInstance(homog)
- 
- Cv = vacancyfluxState(homog)%subState0(1,offset)
- call vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el)
- Cv = Cv + subdt*Cvdot 
- 
- vacancyflux_isochempot_updateState = [     abs(Cv - vacancyfluxState(homog)%state(1,offset)) &
-                                         <= err_vacancyflux_tolAbs &
-                                       .or. abs(Cv - vacancyfluxState(homog)%state(1,offset)) &
-                                         <= err_vacancyflux_tolRel*abs(vacancyfluxState(homog)%state(1,offset)), &
-                                      .true.]
-
- vacancyConc    (homog)%p(vacancyfluxMapping(homog)%p(ip,el)) = Cv
- vacancyConcRate(homog)%p(vacancyfluxMapping(homog)%p(ip,el)) = &
-   (vacancyfluxState(homog)%state(1,offset) - vacancyfluxState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal))
- 
-end function vacancyflux_isochempot_updateState
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  calculates homogenized vacancy driving forces  
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el)
- use material, only: &
-   homogenization_Ngrains, &
-   mappingHomogenization, &
-   phaseAt, &
-   phase_source, &
-   phase_Nsources, &
-   SOURCE_vacancy_phenoplasticity_ID, &
-   SOURCE_vacancy_irradiation_ID, &
-   SOURCE_vacancy_thermalfluc_ID
- use source_vacancy_phenoplasticity, only: &
-   source_vacancy_phenoplasticity_getRateAndItsTangent
- use source_vacancy_irradiation, only: &
-   source_vacancy_irradiation_getRateAndItsTangent
- use source_vacancy_thermalfluc, only: &
-   source_vacancy_thermalfluc_getRateAndItsTangent
- 
- implicit none
- integer(pInt), intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
- real(pReal),   intent(in) :: &
-   Cv
- integer(pInt) :: &
-   phase, &
-   grain, &
-   source
- real(pReal) :: &
-   CvDot, dCvDot_dCv, localCvDot, dLocalCvDot_dCv  
-
- CvDot = 0.0_pReal
- dCvDot_dCv = 0.0_pReal
- do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el))
-   phase = phaseAt(grain,ip,el)
-   do source = 1_pInt, phase_Nsources(phase)
-     select case(phase_source(source,phase))                                                   
-       case (SOURCE_vacancy_phenoplasticity_ID)
-        call source_vacancy_phenoplasticity_getRateAndItsTangent  (localCvDot, dLocalCvDot_dCv, grain, ip,  el)
-
-       case (SOURCE_vacancy_irradiation_ID)
-        call source_vacancy_irradiation_getRateAndItsTangent  (localCvDot, dLocalCvDot_dCv, grain, ip,  el)
-
-       case (SOURCE_vacancy_thermalfluc_ID)
-        call source_vacancy_thermalfluc_getRateAndItsTangent(localCvDot, dLocalCvDot_dCv, grain, ip,  el)
-
-     end select
-     CvDot = CvDot + localCvDot
-     dCvDot_dCv = dCvDot_dCv + dLocalCvDot_dCv
-   enddo  
- enddo
- 
- CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
- dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
- 
-end subroutine vacancyflux_isochempot_getSourceAndItsTangent
-
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of vacancy transport results
-!--------------------------------------------------------------------------------------------------
-function vacancyflux_isochempot_postResults(ip,el)
- use material, only: &
-   mappingHomogenization, &
-   vacancyflux_typeInstance, &
-   vacancyConc, &
-   vacancyfluxMapping
-
- implicit none
- integer(pInt),              intent(in) :: &
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal), dimension(vacancyflux_isochempot_sizePostResults(vacancyflux_typeInstance(mappingHomogenization(2,ip,el)))) :: &
-   vacancyflux_isochempot_postResults
-
- integer(pInt) :: &
-   instance, homog, offset, o, c
-   
- homog     = mappingHomogenization(2,ip,el)
- offset    = vacancyfluxMapping(homog)%p(ip,el)
- instance  = vacancyflux_typeInstance(homog)
-
- c = 0_pInt
- vacancyflux_isochempot_postResults = 0.0_pReal
-
- do o = 1_pInt,vacancyflux_isochempot_Noutput(instance)
-    select case(vacancyflux_isochempot_outputID(o,instance))
- 
-      case (vacancyconc_ID)
-        vacancyflux_isochempot_postResults(c+1_pInt) = vacancyConc(homog)%p(offset)
-        c = c + 1
-    end select
- enddo
-end function vacancyflux_isochempot_postResults
-
-end module vacancyflux_isochempot
--- a/src/vacancyflux_isoconc.f90
+++ b/src/vacancyflux_isoconc.f90
@ -1,62 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for constant vacancy concentration
-!--------------------------------------------------------------------------------------------------
-module vacancyflux_isoconc
-
- implicit none
- private
- 
- public :: &
-   vacancyflux_isoconc_init
-
-contains
-
-!--------------------------------------------------------------------------------------------------
-!> @brief allocates all neccessary fields, reads information from material configuration file
-!--------------------------------------------------------------------------------------------------
-subroutine vacancyflux_isoconc_init()
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
- use, intrinsic :: iso_fortran_env, only: &
-   compiler_version, &
-   compiler_options
-#endif
- use prec, only: &
-   pReal, &
-   pInt 
- use IO, only: &
-   IO_timeStamp
- use material
- use config
- 
- implicit none
- integer(pInt) :: &
-   homog, &
-   NofMyHomog
-
- write(6,'(/,a)')   ' <<<+-  vacancyflux_'//VACANCYFLUX_isoconc_label//' init  -+>>>'
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
-
- initializeInstances: do homog = 1_pInt, material_Nhomogenization
-   
-   myhomog: if (vacancyflux_type(homog) == VACANCYFLUX_isoconc_ID) then
-     NofMyHomog = count(material_homog == homog)
-     vacancyfluxState(homog)%sizeState = 0_pInt
-     vacancyfluxState(homog)%sizePostResults = 0_pInt
-     allocate(vacancyfluxState(homog)%state0   (0_pInt,NofMyHomog))
-     allocate(vacancyfluxState(homog)%subState0(0_pInt,NofMyHomog))
-     allocate(vacancyfluxState(homog)%state    (0_pInt,NofMyHomog))
-     
-     deallocate(vacancyConc    (homog)%p)
-     allocate  (vacancyConc    (homog)%p(1), source=vacancyflux_initialCv(homog))
-     deallocate(vacancyConcRate(homog)%p)
-     allocate  (vacancyConcRate(homog)%p(1), source=0.0_pReal)
-
-   endif myhomog
- enddo initializeInstances
-
-
-end subroutine vacancyflux_isoconc_init
-
-end module vacancyflux_isoconc