diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90
index 3be9032f3..4d8f333ec 100644
--- a/src/Marc/materialpoint_Marc.f90
+++ b/src/Marc/materialpoint_Marc.f90
@@ -158,9 +158,6 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
     else validCalculation
       call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, &
                                                  (elCP-1)*discretization_nIPs + ip)
-      if (.not. terminallyIll) &
-        call homogenization_mechanical_response2(dt,(elCP-1)*discretization_nIPs + ip, &
-                                                    (elCP-1)*discretization_nIPs + ip)
 
       terminalIllness: if (terminallyIll) then
 
diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90
index dae4cd4ee..600c89309 100644
--- a/src/grid/grid_mech_utilities.f90
+++ b/src/grid/grid_mech_utilities.f90
@@ -136,10 +136,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   homogenization_F  = reshape(F,[3,3,product(cells(1:2))*cells3])                                   ! set materialpoint target F to estimated field
 
   call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3)                     ! calculate P field
-  if (.not. terminallyIll) &
-    call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
-  if (.not. terminallyIll) &
-    call homogenization_mechanical_response2(Delta_t,1,product(cells(1:2))*cells3)
 
   P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
   P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index c8e638207..ae672d002 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -323,6 +323,8 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
   real(pREAL), dimension(3,cells(1),cells(2),cells3) :: vectorField
 
 
+  call homogenization_thermal_response(Delta_t_,1,product(cells(1:2))*cells3)
+
   associate(T => x_scal)
     vectorField = utilities_ScalarGradient(T)
     ce = 0
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 726c2fc2e..2da1d73b2 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -168,7 +168,6 @@ module homogenization
   public ::  &
     homogenization_init, &
     homogenization_mechanical_response, &
-    homogenization_mechanical_response2, &
     homogenization_thermal_response, &
     homogenization_thermal_active, &
     homogenization_mu_T, &
@@ -227,7 +226,8 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
     doneAndHappy
 
 
-  !$OMP PARALLEL DO PRIVATE(en,ho,co,converged,doneAndHappy)
+  !$OMP PARALLEL
+  !$OMP DO PRIVATE(en,ho,co,converged,doneAndHappy)
   do ce = cell_start, cell_end
 
     en = material_entry_homogenization(ce)
@@ -260,7 +260,18 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
       terminallyIll = .true.
     end if
   end do
-  !$OMP END PARALLEL DO
+  !$OMP END DO
+
+  !$OMP DO PRIVATE(ho)
+  do ce = cell_start, cell_end
+      ho = material_ID_homogenization(ce)
+      do co = 1, homogenization_Nconstituents(ho)
+        call crystallite_orientations(co,ce)
+      end do
+      call mechanical_homogenize(Delta_t,ce)
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
 
 end subroutine homogenization_mechanical_response
 
@@ -294,32 +305,6 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
 end subroutine homogenization_thermal_response
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief
-!--------------------------------------------------------------------------------------------------
-subroutine homogenization_mechanical_response2(Delta_t,cell_start,cell_end)
-
-  real(pREAL), intent(in) :: Delta_t                                                                !< time increment
-  integer, intent(in) :: &
-    cell_start, cell_end
-
-  integer :: &
-    co, ce, ho
-
-
-  !$OMP PARALLEL DO PRIVATE(ho)
-  do ce = cell_start, cell_end
-      ho = material_ID_homogenization(ce)
-      do co = 1, homogenization_Nconstituents(ho)
-        call crystallite_orientations(co,ce)
-      end do
-      call mechanical_homogenize(Delta_t,ce)
-  end do
-  !$OMP END PARALLEL DO
-
-end subroutine homogenization_mechanical_response2
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes homogenization results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90
index 0764f3443..4b6621b75 100644
--- a/src/mesh/FEM_utilities.f90
+++ b/src/mesh/FEM_utilities.f90
@@ -18,6 +18,7 @@ module FEM_utilities
   use math
   use misc
   use IO
+  use parallelization
   use discretization_mesh
   use homogenization
   use FEM_quadrature
@@ -144,16 +145,15 @@ subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData)
 
   integer(MPI_INTEGER_KIND) :: err_MPI
 
+
   print'(/,1x,a)', '... evaluating constitutive response ......................................'
 
   call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems)                     ! calculate P field
-  if (.not. terminallyIll) &
-    call homogenization_mechanical_response2(Delta_t,1,mesh_maxNips*mesh_NcpElems)
   cutBack = .false.
 
   P_av = sum(homogenization_P,dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+  call parallelization_chkerr(err_MPI)
 
 
 end subroutine utilities_constitutiveResponse