diff --git a/src/constitutive.f90 b/src/constitutive.f90
index e05336cc2..6cf7b5e46 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -120,10 +120,6 @@ module constitutive
       integer, intent(in) :: ph, me
     end subroutine mech_initializeRestorationPoints
 
-    module subroutine constitutive_thermal_initializeRestorationPoints(ph,me)
-      integer, intent(in) :: ph, me
-    end subroutine constitutive_thermal_initializeRestorationPoints
-
 
     module subroutine mech_windForward(ph,me)
       integer, intent(in) :: ph, me
@@ -364,7 +360,6 @@ module constitutive
     constitutive_mech_setF, &
     constitutive_mech_getF, &
     constitutive_initializeRestorationPoints, &
-    constitutive_thermal_initializeRestorationPoints, &
     constitutive_windForward, &
     KINEMATICS_UNDEFINED_ID ,&
     KINEMATICS_CLEAVAGE_OPENING_ID, &
@@ -631,9 +626,6 @@ subroutine crystallite_init()
     do so = 1, phase_Nsources(ph)
       allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0)                 ! ToDo: hack
     enddo
-    do so = 1, thermal_Nsources(ph)
-      allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0)               ! ToDo: hack
-    enddo
   enddo
 
   print'(a42,1x,i10)', '    # of elements:                       ', eMax
diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90
index fe07d2328..56e03ac6f 100644
--- a/src/constitutive_mech.f90
+++ b/src/constitutive_mech.f90
@@ -1586,9 +1586,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   do so = 1, phase_Nsources(ph)
     damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me)
   enddo
-  do so = 1, thermal_Nsources(ph)
-    thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me)
-  enddo
   subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
   subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
   subF0  = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)
diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90
index 610b0b4a6..a82262e81 100644
--- a/src/constitutive_thermal.f90
+++ b/src/constitutive_thermal.f90
@@ -116,8 +116,8 @@ module subroutine thermal_init(phases)
   PhaseLoop2:do ph = 1,phases%length
 
     do so = 1,thermal_Nsources(ph)
-      thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0
-      thermalState(ph)%p(so)%state             = thermalState(ph)%p(so)%partitionedState0
+      deallocate(thermalState(ph)%p(so)%partitionedState0)
+      thermalState(ph)%p(so)%state             = thermalState(ph)%p(so)%state0
     enddo
 
     thermal_source_maxSizeDotState   = max(thermal_source_maxSizeDotState, &
@@ -210,9 +210,6 @@ module function thermal_stress(Delta_t,ph,me) result(converged_)
   integer :: so
 
 
-  do so = 1, thermal_Nsources(ph)
-    thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me)
-  enddo
   converged_ = .not. integrateThermalState(Delta_t,ph,me)
 
 end function thermal_stress
@@ -237,28 +234,13 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
 
   do so = 1, thermal_Nsources(ph)
     sizeDotState = thermalState(ph)%p(so)%sizeDotState
-    thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) &
+    thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%state0(1:sizeDotState,me) &
                                                     + thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t
   enddo
 
 end function integrateThermalState
 
 
-module subroutine constitutive_thermal_initializeRestorationPoints(ph,me)
-
-  integer, intent(in) :: ph, me
-
-  integer :: so
-
-
-  do so = 1, size(thermalState(ph)%p)
-    thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me)
-  enddo
-
-end subroutine constitutive_thermal_initializeRestorationPoints
-
-
-
 module subroutine thermal_forward()
 
   integer :: ph, so
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 627c7c36e..86caafef3 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -279,7 +279,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
         call thermal_partition(homogenization_T(ce),ce)
         do co = 1, homogenization_Nconstituents(ho)
           ph = material_phaseAt(co,el)
-          call constitutive_thermal_initializeRestorationPoints(ph,material_phaseMemberAt(co,ip,el))
           if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
             if (.not. terminallyIll) &                                                           ! so first signals terminally ill...
               print*, ' Integration point ', ip,' at element ', el, ' terminally ill'