From 95c64f7a0aefaebcfa6c82802a6205bba5a19a22 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 21:53:00 +0100
Subject: [PATCH 01/15] avoid failing self test

increase number of samples to have less corner cases.
Needs to be allocatable to avoid stack/heap issue on ifort
---
 src/math.f90 | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/math.f90 b/src/math.f90
index 34cb6eecb..3ea4d6a08 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -1432,9 +1432,11 @@ subroutine selfTest
     error stop 'math_LeviCivita'
 
   normal_distribution: block
-    real(pReal), dimension(500000) :: r
+    integer, parameter :: N = 1000000
+    real(pReal), dimension(:), allocatable :: r
     real(pReal) :: mu, sigma
 
+    allocate(r(N))
     call random_number(mu)
     call random_number(sigma)
 
@@ -1443,11 +1445,11 @@ subroutine selfTest
 
     call math_normal(r,mu,sigma)
 
-    if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) &
+    if (abs(mu -sum(r)/real(N,pReal))>5.0e-2_pReal) &
       error stop 'math_normal(mu)'
 
-    mu = sum(r)/real(size(r),pReal)
-    if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
+    mu = sum(r)/real(N,pReal)
+    if (abs(sigma**2 -1.0_pReal/real(N-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
       error stop 'math_normal(sigma)'
   end block normal_distribution
 

From de9183af4e78d4663c062d2b0174ab2dc2db0910 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 21:55:49 +0100
Subject: [PATCH 02/15] functions without side-effects are 'pure'

basically all 'getter' functions should be pure
---
 src/phase_mechanical.f90         | 6 +++---
 src/phase_mechanical_elastic.f90 | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index 37e6fe7bf..1917d81c9 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -168,17 +168,17 @@ submodule(phase) mechanical
       integer,     intent(in) :: ph,en
     end function plastic_dislotwin_homogenizedC
 
-    module function elastic_C66(ph,en) result(C66)
+    pure module function elastic_C66(ph,en) result(C66)
       real(pReal), dimension(6,6) :: C66
       integer,     intent(in) :: ph, en
     end function elastic_C66
 
-    module function elastic_mu(ph,en) result(mu)
+    pure module function elastic_mu(ph,en) result(mu)
       real(pReal) :: mu
       integer, intent(in) :: ph, en
     end function elastic_mu
 
-    module function elastic_nu(ph,en) result(nu)
+    pure module function elastic_nu(ph,en) result(nu)
       real(pReal) :: nu
       integer, intent(in) :: ph, en
     end function elastic_nu
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index acbbac236..24797a880 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -86,7 +86,7 @@ end subroutine elastic_init
 !--------------------------------------------------------------------------------------------------
 !> @brief return 6x6 elasticity tensor
 !--------------------------------------------------------------------------------------------------
-module function elastic_C66(ph,en) result(C66)
+pure module function elastic_C66(ph,en) result(C66)
 
   integer, intent(in) :: &
     ph, &
@@ -140,7 +140,7 @@ end function elastic_C66
 !--------------------------------------------------------------------------------------------------
 !> @brief return shear modulus
 !--------------------------------------------------------------------------------------------------
-module function elastic_mu(ph,en) result(mu)
+pure module function elastic_mu(ph,en) result(mu)
 
   integer, intent(in) :: &
     ph, &
@@ -157,7 +157,7 @@ end function elastic_mu
 !--------------------------------------------------------------------------------------------------
 !> @brief return Poisson ratio
 !--------------------------------------------------------------------------------------------------
-module function elastic_nu(ph,en) result(nu)
+pure module function elastic_nu(ph,en) result(nu)
 
   integer, intent(in) :: &
     ph, &

From 1661b815b21db384545f5476294ce3b4332c3c43 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 22:02:29 +0100
Subject: [PATCH 03/15] correct calculation of temperature dependent stacking
 fault energy

---
 src/phase_mechanical_plastic_dislotwin.f90 | 13 +++++++------
 src/phase_mechanical_plastic_nonlocal.f90  |  2 +-
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 78681fa24..805b3cc43 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -153,10 +153,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
   print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
-  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'
 
   print'(/,1x,a)', 'F. Roters et al., Computational Materials Science 39:91–95, 2007'
-  print'(  1x,a)', 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.commatsci.2006.04.014'
 
   print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
   print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032'
@@ -656,12 +656,14 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
     dot_gamma_tr
   real(pReal) :: &
     mu, &
-    nu
+    nu, &
+    Gamma
 
   associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
 
     mu = elastic_mu(ph,en)
     nu = elastic_nu(ph,en)
+    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
 
     f_matrix = 1.0_pReal &
              - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
@@ -689,8 +691,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
         else
           ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
           sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
-          b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) &
-                        * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (mu*prm%b_sl(i)), &
+          b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) * Gamma / (mu*prm%b_sl(i)), &
                       1.0_pReal, &
                       prm%ExtendedDislocations)
           v_cl = 2.0_pReal*prm%omega*b_d**2*exp(-prm%Q_cl/(K_B*T)) &
@@ -756,7 +757,7 @@ module subroutine dislotwin_dependentState(T,ph,en)
     sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
     sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T
+    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
 
     !* rescaled volume fraction for topology
     f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw                                              ! this is per system ...
diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90
index 7118055fe..9d88fb5fc 100644
--- a/src/phase_mechanical_plastic_nonlocal.f90
+++ b/src/phase_mechanical_plastic_nonlocal.f90
@@ -203,7 +203,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
   print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
 
   print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
-  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'
 
   print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
   print'(  1x,a)', 'http://publications.rwth-aachen.de/record/229993'

From da00f33487f5020310788273797abd425762713a Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 22:09:53 +0100
Subject: [PATCH 04/15] transformation is only for fcc

---
 src/phase_mechanical_plastic_dislotwin.f90 | 48 ++++++++--------------
 1 file changed, 18 insertions(+), 30 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 805b3cc43..b94a135f1 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -25,9 +25,9 @@ submodule(phase:plastic) dislotwin
       xi_sb               = 1.0_pReal, &                                                            !< value for shearband resistance
       v_sb                = 1.0_pReal, &                                                            !< value for shearband velocity_0
       E_sb                = 1.0_pReal, &                                                            !< activation energy for shear bands
-      delta_G             = 1.0_pReal, &                                                            !< Free energy difference between austensite and martensite
+      delta_G             = 1.0_pReal, &                                                            !< free energy difference between austensite and martensite
       i_tr                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for transformation
-      h                   = 1.0_pReal, &                                                            !< Stack height of hex nucleus
+      h                   = 1.0_pReal, &                                                            !< stack height of hex nucleus
       T_ref               = T_ROOM, &
       a_cI                = 1.0_pReal, &
       a_cF                = 1.0_pReal
@@ -40,14 +40,13 @@ submodule(phase:plastic) dislotwin
       Q_sl,&                                                                                        !< activation energy for glide [J] for each slip system
       v_0, &                                                                                        !< dislocation velocity prefactor [m/s] for each slip system
       dot_N_0_tw, &                                                                                 !< twin nucleation rate [1/m³s] for each twin system
-      dot_N_0_tr, &                                                                                 !< trans nucleation rate [1/m³s] for each trans system
       t_tw, &                                                                                       !< twin thickness [m] for each twin system
       i_sl, &                                                                                       !< Adj. parameter for distance between 2 forest dislocations for each slip system
       t_tr, &                                                                                       !< martensite lamellar thickness [m] for each trans system
       p, &                                                                                          !< p-exponent in glide velocity
       q, &                                                                                          !< q-exponent in glide velocity
-      r, &                                                                                          !< r-exponent in twin nucleation rate
-      s, &                                                                                          !< s-exponent in trans nucleation rate
+      r, &                                                                                          !< exponent in twin nucleation rate
+      s, &                                                                                          !< exponent in trans nucleation rate
       tau_0, &                                                                                      !< strength due to elements in solid solution
       gamma_char, &                                                                                 !< characteristic shear for twins
       B, &                                                                                          !< drag coefficient
@@ -306,10 +305,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%b_tr = pl%get_as1dFloat('b_tr')
       prm%b_tr = math_expand(prm%b_tr,prm%N_tr)
 
-      prm%h             = pl%get_asFloat('h',       defaultVal=0.0_pReal) ! ToDo: How to handle that???
-      prm%i_tr          = pl%get_asFloat('i_tr',    defaultVal=0.0_pReal) ! ToDo: How to handle that???
+      prm%h             = pl%get_asFloat('h',       defaultVal=0.0_pReal) ! ToDo: This is not optional!
+      prm%i_tr          = pl%get_asFloat('i_tr',    defaultVal=0.0_pReal) ! ToDo: This is not optional!
       prm%delta_G       = pl%get_asFloat('delta_G')
-      prm%x_c_tr        = pl%get_asFloat('x_c_tr',  defaultVal=0.0_pReal) ! ToDo: How to handle that???
+      prm%x_c_tr        = pl%get_asFloat('x_c_tr',  defaultVal=0.0_pReal) ! ToDo: This is not optional!
       prm%L_tr          = pl%get_asFloat('L_tr')
       prm%a_cI          = pl%get_asFloat('a_cI',    defaultVal=0.0_pReal)
       prm%a_cF          = pl%get_asFloat('a_cF',    defaultVal=0.0_pReal)
@@ -324,10 +323,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
                                                prm%a_cI, &
                                                prm%a_cF)
 
-      if (phase_lattice(ph) /= 'cF') then
-        prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr')
-        prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr)
-      endif
       prm%t_tr = pl%get_as1dFloat('t_tr')
       prm%t_tr = math_expand(prm%t_tr,prm%N_tr)
       prm%s    = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal])
@@ -339,11 +334,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
       if (    prm%i_tr          < 0.0_pReal)  extmsg = trim(extmsg)//' i_tr'
       if (any(prm%t_tr          < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr'
       if (any(prm%s             < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr'
-      if (phase_lattice(ph) /= 'cF') then
-        if (any(prm%dot_N_0_tr  < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr'
-      end if
     else transActive
-      allocate(prm%s,prm%b_tr,prm%t_tr,prm%dot_N_0_tr,source=emptyRealArray)
+      allocate(prm%s,prm%b_tr,prm%t_tr,source=emptyRealArray)
       allocate(prm%h_tr_tr(0,0))
     end if transActive
 
@@ -1042,20 +1034,16 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
     do i = 1, prm%sum_N_tr
       tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-      isFCC: if (prm%fccTwinTransNucleation) then
-        s1=prm%fcc_twinNucleationSlipPair(1,i)
-        s2=prm%fcc_twinNucleationSlipPair(2,i)
-        if (tau(i) < dst%tau_r_tr(i,en)) then                                                       ! ToDo: correct?
-          Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-                 abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
-                  (prm%L_tr*prm%b_sl(i))*&
-                  (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i))))
-        else
-          Ndot0=0.0_pReal
-        end if
-      else isFCC
-        Ndot0=prm%dot_N_0_tr(i)
-      end if isFCC
+      s1=prm%fcc_twinNucleationSlipPair(1,i)
+      s2=prm%fcc_twinNucleationSlipPair(2,i)
+      if (tau(i) < dst%tau_r_tr(i,en)) then                                                         ! ToDo: correct?
+        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
+               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
+                (prm%L_tr*prm%b_sl(i))*&
+                (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i))))
+      else
+        Ndot0=0.0_pReal
+      end if
     end do
 
     significantStress: where(tau > tol_math_check)

From d130225c9f90939676bb4ff3cfd66420574b2c23 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 22:13:48 +0100
Subject: [PATCH 05/15] polishing

---
 src/phase_mechanical_plastic_dislotwin.f90 | 42 +++++++++++-----------
 src/phase_mechanical_plastic_nonlocal.f90  |  1 -
 2 files changed, 21 insertions(+), 22 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index b94a135f1..09c7ca35c 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -779,9 +779,8 @@ module subroutine dislotwin_dependentState(T,ph,en)
                          + 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) &
                          + prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
 
-    dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2
-    dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2
-
+    dst%V_tw(:,en) = PI/4.0_pReal*dst%Lambda_tw(:,en)**2*prm%t_tw
+    dst%V_tr(:,en) = PI/4.0_pReal*dst%Lambda_tr(:,en)**2*prm%t_tr
 
     x0 = mu*prm%b_tw**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                  ! ToDo: In the paper, this is the Burgers vector for slip
     dst%tau_r_tw(:,en) = mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
@@ -954,11 +953,12 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
   real, dimension(param(ph)%sum_N_tw) :: &
     tau, &
-    Ndot0, &
+    dot_N_0, &
     stressRatio_r, &
     ddot_gamma_dtau
-
-  integer :: i,s1,s2
+  integer, dimension(2) :: &
+    s
+  integer :: i
 
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
@@ -966,24 +966,23 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
     do i = 1, prm%sum_N_tw
       tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
       isFCC: if (prm%fccTwinTransNucleation) then
-        s1=prm%fcc_twinNucleationSlipPair(1,i)
-        s2=prm%fcc_twinNucleationSlipPair(2,i)
         if (tau(i) < dst%tau_r_tw(i,en)) then                                                       ! ToDo: correct?
-          Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-                 abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
+          s=prm%fcc_twinNucleationSlipPair(1:2,i)
+          dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
+                   abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
                   (prm%L_tw*prm%b_sl(i))*&
                   (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tw(i,en)-tau(i))))
         else
-          Ndot0=0.0_pReal
+          dot_N_0=0.0_pReal
         end if
       else isFCC
-        Ndot0=prm%dot_N_0_tw(i)
+        dot_N_0=prm%dot_N_0_tw(i)
       end if isFCC
     end do
 
     significantStress: where(tau > tol_math_check)
       StressRatio_r   = (dst%tau_hat_tw(:,en)/tau)**prm%r
-      dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r)
+      dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,en) * dot_N_0*exp(-StressRatio_r)
       ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
     else where significantStress
       dot_gamma_tw    = 0.0_pReal
@@ -1024,31 +1023,32 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
   real, dimension(param(ph)%sum_N_tr) :: &
     tau, &
-    Ndot0, &
+    dot_N_0, &
     stressRatio_s, &
     ddot_gamma_dtau
-  integer :: i,s1,s2
+  integer, dimension(2) :: &
+    s
+  integer :: i
 
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
     do i = 1, prm%sum_N_tr
       tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-      s1=prm%fcc_twinNucleationSlipPair(1,i)
-      s2=prm%fcc_twinNucleationSlipPair(2,i)
       if (tau(i) < dst%tau_r_tr(i,en)) then                                                         ! ToDo: correct?
-        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
+        s=prm%fcc_twinNucleationSlipPair(1:2,i)
+        dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
+                 abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
                 (prm%L_tr*prm%b_sl(i))*&
                 (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i))))
       else
-        Ndot0=0.0_pReal
+        dot_N_0=0.0_pReal
       end if
     end do
 
     significantStress: where(tau > tol_math_check)
       StressRatio_s   = (dst%tau_hat_tr(:,en)/tau)**prm%s
-      dot_gamma_tr    = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s)
+      dot_gamma_tr    = dst%V_tr(:,en) * dot_N_0*exp(-StressRatio_s)
       ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
     else where significantStress
       dot_gamma_tr  = 0.0_pReal
diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90
index 9d88fb5fc..1715af2d9 100644
--- a/src/phase_mechanical_plastic_nonlocal.f90
+++ b/src/phase_mechanical_plastic_nonlocal.f90
@@ -1570,7 +1570,6 @@ subroutine stateInit(ini,phase,Nentries)
     upto, &
     s
   real(pReal), dimension(2) :: &
-    noise, &
     rnd
   real(pReal) :: &
     meanDensity, &

From f2b6ddece1f44ffb23e587a362cef195e6c890a4 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 22:36:06 +0100
Subject: [PATCH 06/15] reduce memory footprint

---
 src/phase_mechanical_plastic_dislotwin.f90 | 42 ++++++++++++++--------
 1 file changed, 27 insertions(+), 15 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 09c7ca35c..e96c0dbbf 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -103,9 +103,7 @@ submodule(phase:plastic) dislotwin
       tau_hat_tw, &                                                                                 !< threshold stress for twinning
       tau_hat_tr, &                                                                                 !< threshold stress for transformation
       V_tw, &                                                                                       !< volume of a new twin
-      V_tr, &                                                                                       !< volume of a new martensite disc
-      tau_r_tw, &                                                                                   !< stress to bring partials close together (twin)
-      tau_r_tr                                                                                      !< stress to bring partials close together (trans)
+      V_tr                                                                                          !< volume of a new martensite disc
   end type tDislotwinDependentState
 
 !--------------------------------------------------------------------------------------------------
@@ -435,12 +433,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
     allocate(dst%Lambda_tw             (prm%sum_N_tw,Nmembers),source=0.0_pReal)
     allocate(dst%tau_hat_tw            (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_r_tw              (prm%sum_N_tw,Nmembers),source=0.0_pReal)
     allocate(dst%V_tw                  (prm%sum_N_tw,Nmembers),source=0.0_pReal)
 
     allocate(dst%Lambda_tr             (prm%sum_N_tr,Nmembers),source=0.0_pReal)
     allocate(dst%tau_hat_tr            (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_r_tr              (prm%sum_N_tr,Nmembers),source=0.0_pReal)
     allocate(dst%V_tr                  (prm%sum_N_tr,Nmembers),source=0.0_pReal)
 
     end associate
@@ -782,12 +778,6 @@ module subroutine dislotwin_dependentState(T,ph,en)
     dst%V_tw(:,en) = PI/4.0_pReal*dst%Lambda_tw(:,en)**2*prm%t_tw
     dst%V_tr(:,en) = PI/4.0_pReal*dst%Lambda_tr(:,en)**2*prm%t_tr
 
-    x0 = mu*prm%b_tw**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                  ! ToDo: In the paper, this is the Burgers vector for slip
-    dst%tau_r_tw(:,en) = mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
-
-    x0 = mu*prm%b_tr**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                  ! ToDo: In the paper, this is the Burgers vector for slip
-    dst%tau_r_tr(:,en) = mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
-
   end associate
 
 end subroutine dislotwin_dependentState
@@ -956,6 +946,11 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
     dot_N_0, &
     stressRatio_r, &
     ddot_gamma_dtau
+  real :: &
+    x0, &
+    tau_r, &
+    Gamma, &
+    mu, nu
   integer, dimension(2) :: &
     s
   integer :: i
@@ -963,15 +958,21 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
+    mu = elastic_mu(ph,en)
+    nu = elastic_nu(ph,en)
+    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
+
     do i = 1, prm%sum_N_tw
       tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
       isFCC: if (prm%fccTwinTransNucleation) then
-        if (tau(i) < dst%tau_r_tw(i,en)) then                                                       ! ToDo: correct?
+        x0 = mu*prm%b_tw(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                   ! ToDo: In the paper, this is the Burgers vector for slip
+        tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
+        if (tau(i) < tau_r) then                                                                    ! ToDo: correct?
           s=prm%fcc_twinNucleationSlipPair(1:2,i)
           dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
                    abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
                   (prm%L_tw*prm%b_sl(i))*&
-                  (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tw(i,en)-tau(i))))
+                  (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
         else
           dot_N_0=0.0_pReal
         end if
@@ -1026,6 +1027,11 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
     dot_N_0, &
     stressRatio_s, &
     ddot_gamma_dtau
+  real :: &
+    x0, &
+    tau_r, &
+    Gamma, &
+    mu, nu
   integer, dimension(2) :: &
     s
   integer :: i
@@ -1033,14 +1039,20 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
+    mu = elastic_mu(ph,en)
+    nu = elastic_nu(ph,en)
+    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
+
     do i = 1, prm%sum_N_tr
       tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-      if (tau(i) < dst%tau_r_tr(i,en)) then                                                         ! ToDo: correct?
+      x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, this is the Burgers vector for slip              ! ToDo: In the paper, this is the Burgers vector for slip
+      tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
+      if (tau(i) < tau_r) then                                                                      ! ToDo: correct?
         s=prm%fcc_twinNucleationSlipPair(1:2,i)
         dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
                  abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
                 (prm%L_tr*prm%b_sl(i))*&
-                (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i))))
+                (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
       else
         dot_N_0=0.0_pReal
       end if

From 0116e2dae605a062d0f6586a9a18239b7cf208df Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 Dec 2021 22:38:39 +0100
Subject: [PATCH 07/15] bugfix: write only to active twin/trans system

---
 src/phase_mechanical_plastic_dislotwin.f90 | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index e96c0dbbf..b8f3516ed 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -969,15 +969,14 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
         tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
         if (tau(i) < tau_r) then                                                                    ! ToDo: correct?
           s=prm%fcc_twinNucleationSlipPair(1:2,i)
-          dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                   abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
-                  (prm%L_tw*prm%b_sl(i))*&
-                  (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
+          dot_N_0(i)=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
+                      abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
+                     (prm%L_tw*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
         else
-          dot_N_0=0.0_pReal
+          dot_N_0(i)=0.0_pReal
         end if
       else isFCC
-        dot_N_0=prm%dot_N_0_tw(i)
+        dot_N_0(i)=prm%dot_N_0_tw(i)
       end if isFCC
     end do
 
@@ -1045,16 +1044,15 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
     do i = 1, prm%sum_N_tr
       tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-      x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, this is the Burgers vector for slip              ! ToDo: In the paper, this is the Burgers vector for slip
+      x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, this is the Burgers vector for slip
       tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
       if (tau(i) < tau_r) then                                                                      ! ToDo: correct?
         s=prm%fcc_twinNucleationSlipPair(1:2,i)
-        dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                 abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
-                (prm%L_tr*prm%b_sl(i))*&
-                (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
+        dot_N_0(i)=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
+                    abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
+                   (prm%L_tr*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
       else
-        dot_N_0=0.0_pReal
+        dot_N_0(i)=0.0_pReal
       end if
     end do
 

From 017c182640a079b2b4cfdfc37d35506d95d0a0f7 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 31 Dec 2021 08:45:55 +0100
Subject: [PATCH 08/15] branch only once

---
 src/phase_mechanical_plastic_dislotwin.f90 | 36 +++++++++-------------
 1 file changed, 15 insertions(+), 21 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index b8f3516ed..5e4b757c7 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -966,7 +966,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
       tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
       isFCC: if (prm%fccTwinTransNucleation) then
         x0 = mu*prm%b_tw(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                   ! ToDo: In the paper, this is the Burgers vector for slip
-        tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
+        tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0)
         if (tau(i) < tau_r) then                                                                    ! ToDo: correct?
           s=prm%fcc_twinNucleationSlipPair(1:2,i)
           dot_N_0(i)=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
@@ -1022,13 +1022,12 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
     ddot_gamma_dtau_tr
 
   real, dimension(param(ph)%sum_N_tr) :: &
-    tau, &
-    dot_N_0, &
-    stressRatio_s, &
     ddot_gamma_dtau
   real :: &
+    stressRatio_s, &
+    tau, tau_r, &
+    dot_N_0, &
     x0, &
-    tau_r, &
     Gamma, &
     mu, nu
   integer, dimension(2) :: &
@@ -1043,28 +1042,23 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
     Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
 
     do i = 1, prm%sum_N_tr
-      tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
+      tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
       x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, this is the Burgers vector for slip
-      tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
-      if (tau(i) < tau_r) then                                                                      ! ToDo: correct?
+      tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0)
+      if (tau > tol_math_check .and. tau < tau_r) then
         s=prm%fcc_twinNucleationSlipPair(1:2,i)
-        dot_N_0(i)=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                    abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
-                   (prm%L_tr*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
+        dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
+                 abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
+                (prm%L_tr*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)))
+        StressRatio_s = (dst%tau_hat_tr(i,en)/tau)**prm%s(i)
+        dot_gamma_tr(i) = dst%V_tr(i,en) * dot_N_0*exp(-StressRatio_s)
+        ddot_gamma_dtau(i) = (dot_gamma_tr(i)*prm%s(i)/tau)*StressRatio_s
       else
-        dot_N_0(i)=0.0_pReal
+        dot_gamma_tr(i) = 0.0_pReal
+        ddot_gamma_dtau(i) = 0.0_pReal
       end if
     end do
 
-    significantStress: where(tau > tol_math_check)
-      StressRatio_s   = (dst%tau_hat_tr(:,en)/tau)**prm%s
-      dot_gamma_tr    = dst%V_tr(:,en) * dot_N_0*exp(-StressRatio_s)
-      ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
-    else where significantStress
-      dot_gamma_tr  = 0.0_pReal
-      ddot_gamma_dtau = 0.0_pReal
-    end where significantStress
-
   end associate
 
   if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau

From e678b231d99859f098e12f65f97b5214b2ee3608 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 31 Dec 2021 09:01:23 +0100
Subject: [PATCH 09/15] following naming convention

---
 src/phase_mechanical_plastic_dislotwin.f90 | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 5e4b757c7..3c666934a 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -944,7 +944,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
   real, dimension(param(ph)%sum_N_tw) :: &
     tau, &
     dot_N_0, &
-    stressRatio_r, &
+    ratio_tau_r, &
     ddot_gamma_dtau
   real :: &
     x0, &
@@ -981,11 +981,11 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
     end do
 
     significantStress: where(tau > tol_math_check)
-      StressRatio_r   = (dst%tau_hat_tw(:,en)/tau)**prm%r
-      dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,en) * dot_N_0*exp(-StressRatio_r)
-      ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
+      ratio_tau_r = (dst%tau_hat_tw(:,en)/tau)**prm%r
+      dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * dot_N_0*exp(-ratio_tau_r)
+      ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*ratio_tau_r
     else where significantStress
-      dot_gamma_tw    = 0.0_pReal
+      dot_gamma_tw = 0.0_pReal
       ddot_gamma_dtau = 0.0_pReal
     end where significantStress
 
@@ -1024,7 +1024,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
   real, dimension(param(ph)%sum_N_tr) :: &
     ddot_gamma_dtau
   real :: &
-    stressRatio_s, &
+    ratio_tau_s, &
     tau, tau_r, &
     dot_N_0, &
     x0, &
@@ -1046,13 +1046,16 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
       x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, this is the Burgers vector for slip
       tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0)
       if (tau > tol_math_check .and. tau < tau_r) then
+
+        ratio_tau_s = (dst%tau_hat_tr(i,en)/tau)**prm%s(i)
+
         s=prm%fcc_twinNucleationSlipPair(1:2,i)
         dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
                  abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
                 (prm%L_tr*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)))
-        StressRatio_s = (dst%tau_hat_tr(i,en)/tau)**prm%s(i)
-        dot_gamma_tr(i) = dst%V_tr(i,en) * dot_N_0*exp(-StressRatio_s)
-        ddot_gamma_dtau(i) = (dot_gamma_tr(i)*prm%s(i)/tau)*StressRatio_s
+
+        dot_gamma_tr(i) = dst%V_tr(i,en) * dot_N_0*exp(-ratio_tau_s)
+        ddot_gamma_dtau(i) = (dot_gamma_tr(i)*prm%s(i)/tau)*ratio_tau_s
       else
         dot_gamma_tr(i) = 0.0_pReal
         ddot_gamma_dtau(i) = 0.0_pReal

From 770cf33667739b80daa68b170f666d2f16d9b66f Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 1 Jan 2022 11:37:04 +0100
Subject: [PATCH 10/15] correct calculation of tangent.

thanks to Seyedamirhossein Motaman (RWTH Aachen) for reporting
---
 src/phase_mechanical_plastic_dislotwin.f90 | 116 ++++++++++++---------
 1 file changed, 64 insertions(+), 52 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 3c666934a..5d2191726 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -17,8 +17,8 @@ submodule(phase:plastic) dislotwin
       p_sb                = 1.0_pReal, &                                                            !< p-exponent in shear band velocity
       q_sb                = 1.0_pReal, &                                                            !< q-exponent in shear band velocity
       i_tw                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for twinning
-      L_tw                = 1.0_pReal, &                                                            !< Length of twin nuclei in Burgers vectors
-      L_tr                = 1.0_pReal, &                                                            !< Length of trans nuclei in Burgers vectors
+      L_tw                = 1.0_pReal, &                                                            !< length of twin nuclei in Burgers vectors: TODO unit should be meters
+      L_tr                = 1.0_pReal, &                                                            !< length of trans nuclei in Burgers vectors: TODO unit should be meters
       x_c_tw              = 1.0_pReal, &                                                            !< critical distance for formation of twin nucleus
       x_c_tr              = 1.0_pReal, &                                                            !< critical distance for formation of trans nucleus
       V_cs                = 1.0_pReal, &                                                            !< cross slip volume
@@ -731,8 +731,6 @@ module subroutine dislotwin_dependentState(T,ph,en)
    real(pReal), dimension(param(ph)%sum_N_tr) :: &
     inv_lambda_tr_tr, &                                                                             !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite
     f_over_t_tr
-  real(pReal), dimension(:), allocatable :: &
-    x0
   real(pReal) :: &
     mu, &
     nu
@@ -941,16 +939,15 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
   real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
     ddot_gamma_dtau_tw
 
-  real, dimension(param(ph)%sum_N_tw) :: &
-    tau, &
-    dot_N_0, &
-    ratio_tau_r, &
-    ddot_gamma_dtau
   real :: &
+    ratio_tau_s, &
+    tau, tau_r, &
+    dot_N_0, &
     x0, &
-    tau_r, &
     Gamma, &
-    mu, nu
+    mu, nu, &
+    P_ncs, dP_ncs_dtau, &
+    P, dP_dtau
   integer, dimension(2) :: &
     s
   integer :: i
@@ -958,41 +955,53 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    mu = elastic_mu(ph,en)
-    nu = elastic_nu(ph,en)
-    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
+    isFCC: if (prm%fccTwinTransNucleation) then
+      mu = elastic_mu(ph,en)
+      nu = elastic_nu(ph,en)
+      Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
 
-    do i = 1, prm%sum_N_tw
-      tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
-      isFCC: if (prm%fccTwinTransNucleation) then
-        x0 = mu*prm%b_tw(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                   ! ToDo: In the paper, this is the Burgers vector for slip
+      do i = 1, prm%sum_N_tw
+        tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
+        x0 = mu*prm%b_tw(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                   ! ToDo: In the paper, the Burgers vector for slip is used
         tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0)
-        if (tau(i) < tau_r) then                                                                    ! ToDo: correct?
-          s=prm%fcc_twinNucleationSlipPair(1:2,i)
-          dot_N_0(i)=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                      abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
-                     (prm%L_tw*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau(i))))
-        else
-          dot_N_0(i)=0.0_pReal
-        end if
-      else isFCC
-        dot_N_0(i)=prm%dot_N_0_tw(i)
-      end if isFCC
-    end do
 
-    significantStress: where(tau > tol_math_check)
-      ratio_tau_r = (dst%tau_hat_tw(:,en)/tau)**prm%r
-      dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * dot_N_0*exp(-ratio_tau_r)
-      ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*ratio_tau_r
-    else where significantStress
-      dot_gamma_tw = 0.0_pReal
-      ddot_gamma_dtau = 0.0_pReal
-    end where significantStress
+        if (tau > tol_math_check .and. tau < tau_r) then
+          ratio_tau_s = (dst%tau_hat_tw(i,en)/tau)**prm%r(i)
+          P = exp(-ratio_tau_s)
+          dP_dTau = prm%r(i) * ratio_tau_s/tau * P
+
+          s=prm%fcc_twinNucleationSlipPair(1:2,i)
+          dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
+                   abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/(prm%L_tw*prm%b_sl(i))
+
+          P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
+          dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
+
+          dot_gamma_tw(i) = dst%V_tw(i,en)*dot_N_0*P_ncs*P
+          if (present(ddot_gamma_dtau_tw)) &
+            ddot_gamma_dtau_tw(i) = dst%V_tw(i,en)*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
+        else
+          dot_gamma_tw(i) = 0.0_pReal
+          if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
+        end if
+      end do
+
+    else isFCC
+      do i = 1, prm%sum_N_tw
+        error stop 'not implemented'
+        tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
+        if (tau > tol_math_check) then
+          dot_gamma_tw(i) = 0.0_pReal
+          if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
+        else
+          dot_gamma_tw(i) = 0.0_pReal
+          if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
+        end if
+      end do
+    end if isFCC
 
   end associate
 
-  if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
-
 end subroutine kinetics_tw
 
 
@@ -1021,15 +1030,15 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
   real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
     ddot_gamma_dtau_tr
 
-  real, dimension(param(ph)%sum_N_tr) :: &
-    ddot_gamma_dtau
   real :: &
     ratio_tau_s, &
     tau, tau_r, &
     dot_N_0, &
     x0, &
     Gamma, &
-    mu, nu
+    mu, nu, &
+    P_ncs, dP_ncs_dtau, &
+    P, dP_dtau
   integer, dimension(2) :: &
     s
   integer :: i
@@ -1043,29 +1052,32 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
     do i = 1, prm%sum_N_tr
       tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-      x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, this is the Burgers vector for slip
+      x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, the Burgers vector for slip is used
       tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0)
-      if (tau > tol_math_check .and. tau < tau_r) then
 
+      if (tau > tol_math_check .and. tau < tau_r) then
         ratio_tau_s = (dst%tau_hat_tr(i,en)/tau)**prm%s(i)
+        P = exp(-ratio_tau_s)
+        dP_dTau = prm%s(i) * ratio_tau_s/tau * P
 
         s=prm%fcc_twinNucleationSlipPair(1:2,i)
         dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                 abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/&
-                (prm%L_tr*prm%b_sl(i))*(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)))
+                 abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/(prm%L_tr*prm%b_sl(i))
 
-        dot_gamma_tr(i) = dst%V_tr(i,en) * dot_N_0*exp(-ratio_tau_s)
-        ddot_gamma_dtau(i) = (dot_gamma_tr(i)*prm%s(i)/tau)*ratio_tau_s
+        P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
+        dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
+
+        dot_gamma_tr(i) = dst%V_tr(i,en)*dot_N_0*P_ncs*P
+        if (present(ddot_gamma_dtau_tr)) &
+          ddot_gamma_dtau_tr(i) = dst%V_tr(i,en)*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
       else
         dot_gamma_tr(i) = 0.0_pReal
-        ddot_gamma_dtau(i) = 0.0_pReal
+        if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal
       end if
     end do
 
   end associate
 
-  if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
-
 end subroutine kinetics_tr
 
 end submodule dislotwin

From d181b988c14ba48469658e115424c0b5bc62970d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 1 Jan 2022 10:47:15 +0100
Subject: [PATCH 11/15] using vector access

---
 src/phase_mechanical_plastic_dislotwin.f90 | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 5d2191726..a2f0ce732 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -970,9 +970,9 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
           P = exp(-ratio_tau_s)
           dP_dTau = prm%r(i) * ratio_tau_s/tau * P
 
-          s=prm%fcc_twinNucleationSlipPair(1:2,i)
-          dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                   abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/(prm%L_tw*prm%b_sl(i))
+          s = prm%fcc_twinNucleationSlipPair(1:2,i)
+          dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) &
+                  / (prm%L_tw*prm%b_sl(i))
 
           P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
           dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
@@ -1060,9 +1060,9 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
         P = exp(-ratio_tau_s)
         dP_dTau = prm%s(i) * ratio_tau_s/tau * P
 
-        s=prm%fcc_twinNucleationSlipPair(1:2,i)
-        dot_N_0=(abs(dot_gamma_sl(s(1)))*(stt%rho_mob(s(2),en)+stt%rho_dip(s(2),en))+&
-                 abs(dot_gamma_sl(s(2)))*(stt%rho_mob(s(1),en)+stt%rho_dip(s(1),en)))/(prm%L_tr*prm%b_sl(i))
+        s = prm%fcc_twinNucleationSlipPair(1:2,i)
+        dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) &
+                / (prm%L_tr*prm%b_sl(i))
 
         P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
         dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)

From 03c6708629ce603ff560555aaaa0c5389b1456bd Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 1 Jan 2022 11:26:31 +0100
Subject: [PATCH 12/15] polishing

---
 src/phase_mechanical_plastic_dislotwin.f90 | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index a2f0ce732..6a9043cde 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -940,7 +940,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
     ddot_gamma_dtau_tw
 
   real :: &
-    ratio_tau_s, &
+    ratio_tau_r, &
     tau, tau_r, &
     dot_N_0, &
     x0, &
@@ -962,13 +962,13 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
       do i = 1, prm%sum_N_tw
         tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
-        x0 = mu*prm%b_tw(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                   ! ToDo: In the paper, the Burgers vector for slip is used
-        tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0)
+        x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu))                   ! ToDo: In the paper, the Burgers vector for slip is used
+        tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0)      ! ToDo: In the paper, the Burgers vector for slip is used
 
         if (tau > tol_math_check .and. tau < tau_r) then
-          ratio_tau_s = (dst%tau_hat_tw(i,en)/tau)**prm%r(i)
-          P = exp(-ratio_tau_s)
-          dP_dTau = prm%r(i) * ratio_tau_s/tau * P
+          ratio_tau_r = (dst%tau_hat_tw(i,en)/tau)**prm%r(i)
+          P = exp(-ratio_tau_r)
+          dP_dTau = prm%r(i) * ratio_tau_r/tau * P
 
           s = prm%fcc_twinNucleationSlipPair(1:2,i)
           dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) &
@@ -1052,8 +1052,8 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
     do i = 1, prm%sum_N_tr
       tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-      x0 = mu*prm%b_tr(i)**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                     ! ToDo: In the paper, the Burgers vector for slip is used
-      tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0)
+      x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu))                     ! ToDo: In the paper, the Burgers vector for slip is used
+      tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0)        ! ToDo: In the paper, the Burgers vector for slip is used
 
       if (tau > tol_math_check .and. tau < tau_r) then
         ratio_tau_s = (dst%tau_hat_tr(i,en)/tau)**prm%s(i)

From 4a7f23069c87fe0c4a04970ab36d92312e6dfeff Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 4 Jan 2022 07:32:56 +0100
Subject: [PATCH 13/15] avoid misleading variable name

---
 src/phase_mechanical_plastic_dislotwin.f90 | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 6a9043cde..e1ccf3f5f 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -940,7 +940,6 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
     ddot_gamma_dtau_tw
 
   real :: &
-    ratio_tau_r, &
     tau, tau_r, &
     dot_N_0, &
     x0, &
@@ -966,9 +965,8 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
         tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0)      ! ToDo: In the paper, the Burgers vector for slip is used
 
         if (tau > tol_math_check .and. tau < tau_r) then
-          ratio_tau_r = (dst%tau_hat_tw(i,en)/tau)**prm%r(i)
-          P = exp(-ratio_tau_r)
-          dP_dTau = prm%r(i) * ratio_tau_r/tau * P
+          P = exp(-(dst%tau_hat_tw(i,en)/tau)**prm%r(i))
+          dP_dTau = prm%r(i) * (dst%tau_hat_tw(i,en)/tau)**prm%r(i)/tau * P
 
           s = prm%fcc_twinNucleationSlipPair(1:2,i)
           dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) &
@@ -1031,7 +1029,6 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
     ddot_gamma_dtau_tr
 
   real :: &
-    ratio_tau_s, &
     tau, tau_r, &
     dot_N_0, &
     x0, &
@@ -1056,9 +1053,8 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
       tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0)        ! ToDo: In the paper, the Burgers vector for slip is used
 
       if (tau > tol_math_check .and. tau < tau_r) then
-        ratio_tau_s = (dst%tau_hat_tr(i,en)/tau)**prm%s(i)
-        P = exp(-ratio_tau_s)
-        dP_dTau = prm%s(i) * ratio_tau_s/tau * P
+        P = exp(-(dst%tau_hat_tr(i,en)/tau)**prm%s(i))
+        dP_dTau = prm%s(i) * (dst%tau_hat_tr(i,en)/tau)**prm%s(i)/tau * P
 
         s = prm%fcc_twinNucleationSlipPair(1:2,i)
         dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) &

From 03e3fbd98f0b55d492c434a33cf0decde7f511c4 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 4 Jan 2022 07:40:58 +0100
Subject: [PATCH 14/15] compute only when needed

---
 src/phase_mechanical_plastic_dislotwin.f90 | 25 +++++++++-------------
 1 file changed, 10 insertions(+), 15 deletions(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index e1ccf3f5f..41b4854a7 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -101,9 +101,7 @@ submodule(phase:plastic) dislotwin
       Lambda_tr, &                                                                                  !< mean free path between 2 obstacles seen by a growing martensite
       tau_pass, &                                                                                   !< threshold stress for slip
       tau_hat_tw, &                                                                                 !< threshold stress for twinning
-      tau_hat_tr, &                                                                                 !< threshold stress for transformation
-      V_tw, &                                                                                       !< volume of a new twin
-      V_tr                                                                                          !< volume of a new martensite disc
+      tau_hat_tr                                                                                    !< threshold stress for transformation
   end type tDislotwinDependentState
 
 !--------------------------------------------------------------------------------------------------
@@ -433,11 +431,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
     allocate(dst%Lambda_tw             (prm%sum_N_tw,Nmembers),source=0.0_pReal)
     allocate(dst%tau_hat_tw            (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%V_tw                  (prm%sum_N_tw,Nmembers),source=0.0_pReal)
 
     allocate(dst%Lambda_tr             (prm%sum_N_tr,Nmembers),source=0.0_pReal)
     allocate(dst%tau_hat_tr            (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%V_tr                  (prm%sum_N_tr,Nmembers),source=0.0_pReal)
 
     end associate
 
@@ -773,9 +769,6 @@ module subroutine dislotwin_dependentState(T,ph,en)
                          + 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) &
                          + prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
 
-    dst%V_tw(:,en) = PI/4.0_pReal*dst%Lambda_tw(:,en)**2*prm%t_tw
-    dst%V_tr(:,en) = PI/4.0_pReal*dst%Lambda_tr(:,en)**2*prm%t_tr
-
   end associate
 
 end subroutine dislotwin_dependentState
@@ -942,7 +935,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
   real :: &
     tau, tau_r, &
     dot_N_0, &
-    x0, &
+    x0, V, &
     Gamma, &
     mu, nu, &
     P_ncs, dP_ncs_dtau, &
@@ -975,9 +968,10 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
           P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
           dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
 
-          dot_gamma_tw(i) = dst%V_tw(i,en)*dot_N_0*P_ncs*P
+          V = PI/4.0_pReal*dst%Lambda_tw(i,en)**2*prm%t_tw(i)
+          dot_gamma_tw(i) = V*dot_N_0*P_ncs*P
           if (present(ddot_gamma_dtau_tw)) &
-            ddot_gamma_dtau_tw(i) = dst%V_tw(i,en)*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
+            ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
         else
           dot_gamma_tw(i) = 0.0_pReal
           if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
@@ -1031,7 +1025,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
   real :: &
     tau, tau_r, &
     dot_N_0, &
-    x0, &
+    x0, V, &
     Gamma, &
     mu, nu, &
     P_ncs, dP_ncs_dtau, &
@@ -1062,10 +1056,11 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
         P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
         dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
-
-        dot_gamma_tr(i) = dst%V_tr(i,en)*dot_N_0*P_ncs*P
+        
+        V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i)
+        dot_gamma_tr(i) = V*dot_N_0*P_ncs*P
         if (present(ddot_gamma_dtau_tr)) &
-          ddot_gamma_dtau_tr(i) = dst%V_tr(i,en)*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
+          ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
       else
         dot_gamma_tr(i) = 0.0_pReal
         if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal

From 1140625b1272754102adc6f8a276563cb241eb32 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 4 Jan 2022 17:29:19 +0100
Subject: [PATCH 15/15] copy and paste error

---
 src/phase_mechanical_plastic_dislotwin.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 41b4854a7..bc6fb0ee6 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -954,7 +954,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
       do i = 1, prm%sum_N_tw
         tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
-        x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu))                   ! ToDo: In the paper, the Burgers vector for slip is used
+        x0 = mu*prm%b_tw(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu))                   ! ToDo: In the paper, the Burgers vector for slip is used
         tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0)      ! ToDo: In the paper, the Burgers vector for slip is used
 
         if (tau > tol_math_check .and. tau < tau_r) then