diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90
index c38b9733d..cb6cf0f79 100644
--- a/src/plastic_phenopowerlaw.f90
+++ b/src/plastic_phenopowerlaw.f90
@@ -120,8 +120,6 @@ subroutine plastic_phenopowerlaw_init
    debug_constitutive,&
    debug_levelBasic
  use math, only: &
-   math_Mandel3333to66, &
-   math_Voigt66to3333, &
    math_expand
  use IO, only: &
    IO_warning, &
@@ -487,8 +485,7 @@ pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar99,Mstar,ipc,
    dNeq0
  use math, only: &
    math_mul33xx33,&
-   math_Plain3333to99, &
-   math_Mandel6to33
+   math_Plain3333to99
  use lattice, only: &
    lattice_Sslip, &
    lattice_Stwin, &
@@ -748,6 +745,9 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el)
    plasticState, &
    phasememberAt, &
    phase_plasticityInstance
+ use math, only: &
+   math_mul33xx33, &
+   math_Mandel6to33 
  use lattice, only: &
    lattice_Sslip_v, &
    lattice_Stwin_v, &
@@ -763,6 +763,8 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el)
    ip, &                                                                                            !< integration point
    el                                                                                               !< element                                                                                        !< microstructure state
 
+ real(pReal), dimension(3,3) :: &
+   Mstar
  real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: &
    plastic_phenopowerlaw_postResults
 
@@ -782,6 +784,7 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el)
  associate( prm => param(phase_plasticityInstance(ph)), &
             stt => state(phase_plasticityInstance(ph)), &
             dst => dotState(phase_plasticityInstance(ph)))
+ Mstar = math_Mandel6to33(Mstar6)
 
  plastic_phenopowerlaw_postResults = 0.0_pReal
  c = 0_pInt
@@ -798,23 +801,19 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el)
 
      case (shearrate_slip_ID)
        j = 0_pInt
-       slipFamilies1: do f = 1_pInt,size(prm%Nslip,1)
-         index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph))                                ! at which index starts my family
+       slipFamilies1: do f =1_pInt,size(prm%Nslip,1)
+         index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph))                                   ! at which index starts my family
          slipSystems1: do i = 1_pInt,prm%Nslip(f)
-           j = j + 1_pInt
-           tau_slip_pos  = dot_product(Mstar6,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
+           j = j+1_pInt
+           tau_slip_pos  = math_mul33xx33(Mstar,prm%Schmid_slip(1:3,1:3,j))
            tau_slip_neg  = tau_slip_pos
-           do k = 1,lattice_NnonSchmid(ph)
-             tau_slip_pos = tau_slip_pos +prm%nonSchmidCoeff(k)* &
-                                   dot_product(Mstar6,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph))
-             tau_slip_neg = tau_slip_neg +prm%nonSchmidCoeff(k)* &
-                                   dot_product(Mstar6,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
-           enddo
+           nonSchmidSystems: do k = 1,size(prm%nonSchmidCoeff)
+             tau_slip_pos = tau_slip_pos + math_mul33xx33(Mstar,prm%nonSchmid_pos(1:3,1:3,k,j))
+             tau_slip_neg = tau_slip_neg + math_mul33xx33(Mstar,prm%nonSchmid_neg(1:3,1:3,k,j))
+           enddo nonSchmidSystems
            plastic_phenopowerlaw_postResults(c+j) = prm%gdot0_slip*0.5_pReal* &
-                    ((abs(tau_slip_pos)/stt%s_slip(j,of))**prm%n_slip &
-                    *sign(1.0_pReal,tau_slip_pos) &
-                    +(abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip &
-                    *sign(1.0_pReal,tau_slip_neg))
+                        ( (abs(tau_slip_pos)/(stt%s_slip(j,of)))**prm%n_slip*sign(1.0_pReal,tau_slip_pos) &
+                         +(abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip*sign(1.0_pReal,tau_slip_neg))
          enddo slipSystems1
        enddo slipFamilies1
        c = c + prm%totalNslip
@@ -825,8 +824,7 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el)
          index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph))                                ! at which index starts my family
          slipSystems2: do i = 1_pInt,prm%Nslip(f)
            j = j + 1_pInt
-           plastic_phenopowerlaw_postResults(c+j) = &
-                             dot_product(Mstar6,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
+           plastic_phenopowerlaw_postResults(c+j) = math_mul33xx33(Mstar,prm%Schmid_slip(1:3,1:3,j))
          enddo slipSystems2
        enddo slipFamilies2
        c = c + prm%totalNslip