diff --git a/src/lattice.f90 b/src/lattice.f90
index 090dad385..46f243561 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -1229,6 +1229,7 @@ real(pReal), dimension(4,36), parameter, private :: &
   LATTICE_hex_ID, &
   lattice_SchmidMatrix_slip, &
   lattice_SchmidMatrix_twin, &
+  lattice_nonSchmidMatrix, &
   lattice_interactionSlipSlip2, &
   lattice_interactionTwinTwin2, &
   lattice_interactionSlipTwin2, &
@@ -2226,13 +2227,14 @@ function lattice_C66_trans(Ntrans,C66_parent,structure_parent,cOverA_parent, &
 
 !       ! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction)
 !       ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 (corresponds to their "n1" for positive and negative slip direction respectively)
+
 function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
  use math, only: &
-  INRAD, &
-  math_tensorproduct33, &
-  math_crossproduct, &
-  math_mul33x3, &
-  math_axisAngleToR
+   INRAD, &
+   math_tensorproduct33, &
+   math_crossproduct, &
+   math_mul33x3, &
+   math_axisAngleToR
  implicit none
  integer(pInt),    dimension(:),                   intent(in) :: Nslip                              !< number of active slip systems per family
  real(pReal),      dimension(6),                   intent(in) :: nonSchmidCoefficients
@@ -2253,7 +2255,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
  do i = 1_pInt,sum(Nslip)
    direction = coordinateSystem(1:3,1,i)
    normal    = coordinateSystem(1:3,2,i)
-   np = math_mul33x3(math_axisAngleToR(+direction,60.0_pReal*INRAD), normal)
+   np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal)
      nonSchmidMatrix(1:3,1:3,i) &
        = nonSchmidMatrix(1:3,1:3,i) &
        + nonSchmidCoefficients(1) * math_tensorproduct33(direction, np) &
diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90
index 792889c13..a53f0103f 100644
--- a/src/plastic_phenopowerlaw.f90
+++ b/src/plastic_phenopowerlaw.f90
@@ -213,7 +213,8 @@ subroutine plastic_phenopowerlaw_init
                                                                          defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
      prm%nonSchmidCoeff       = config_phase(p)%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)
      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')
@@ -381,32 +382,9 @@ subroutine plastic_phenopowerlaw_init
      allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal)
 
 
-!--------------------------------------------------------------------------------------------------
-! calculate hardening matrices
-   allocate(prm%nonSchmid_pos(3,3,prm%totalNslip),source = 0.0_pReal)
-   allocate(prm%nonSchmid_neg(3,3,prm%totalNslip),source = 0.0_pReal)
-   i = 0_pInt
-   mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1)                                    ! >>> interaction slip -- X
-     index_myFamily = sum(prm%Nslip(1:f-1_pInt))
-
-     mySlipSystems: do j = 1_pInt,prm%Nslip(f)
-       i = i + 1_pInt
-       prm%nonSchmid_pos(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1,  index_myFamily+j,p)
-       prm%nonSchmid_neg(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1,  index_myFamily+j,p)
-       do k = 1,size(prm%nonSchmidCoeff)
-         prm%nonSchmid_pos(1:3,1:3,i) = prm%nonSchmid_pos(1:3,1:3,i) &
-                                      + lattice_Sslip(1:3,1:3,2*k,  index_myFamily+j,p) &
-                                                                            * prm%nonSchmidCoeff(k)
-         prm%nonSchmid_neg(1:3,1:3,i) = prm%nonSchmid_neg(1:3,1:3,i) &
-                                      + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) &
-                                                                            * prm%nonSchmidCoeff(k)
-       enddo
-     enddo mySlipSystems
-   enddo mySlipFamilies
-
    allocate(prm%gamma_twin_char(prm%totalNtwin),source  = 0.0_pReal)
    i = 0_pInt
-   myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1)                                    ! >>> interaction twin -- X
+   myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1)
      index_myFamily = sum(prm%Ntwin(1:f-1_pInt))
      myTwinSystems: do j = 1_pInt,prm%Ntwin(f)
        i = i + 1_pInt