From 38755cb30d5bcc837235bbb891342e342495b5da Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 16 Aug 2019 13:56:52 -0700
Subject: [PATCH 1/2] polishing

---
 src/plastic_dislotwin.f90 | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90
index e2a8b1728..f1c42b678 100644
--- a/src/plastic_dislotwin.f90
+++ b/src/plastic_dislotwin.f90
@@ -59,7 +59,7 @@ module plastic_dislotwin
      q_sb, &                                                                                        !< q-exponent in shear band velocity
      CEdgeDipMinDistance, &                                                                         !<
      i_tw, &                                                                                        !<
-     tau_0, &                                                                                       !<strength due to elements in solid solution
+     tau_0, &                                                                                       !< strength due to elements in solid solution
      L_tw, &                                                                                        !< Length of twin nuclei in Burgers vectors
      L_tr, &                                                                                        !< Length of trans nuclei in Burgers vectors
      xc_twin, &                                                                                     !< critical distance for formation of twin nucleus
@@ -276,19 +276,16 @@ subroutine plastic_dislotwin_init
      prm%CEdgeDipMinDistance  = config%getFloat('cedgedipmindistance')
      prm%D0                   = config%getFloat('d0')
      prm%Qsd                  = config%getFloat('qsd')
-     prm%omega                = config%getFloat('omega',defaultVal=1000.0_pReal)
      prm%atomicVolume         = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal
-     prm%SFE_0K               = config%getFloat('sfe_0k',defaultVal = 0.0_pReal)
+     prm%SFE_0K               = config%getFloat('sfe_0k', defaultVal = 0.0_pReal)
      prm%dSFE_dT              = config%getFloat('dsfe_dt',defaultVal = 0.0_pReal)
      
      prm%climbSwitch          = config%getInt('climbswitch',defaultVal = 1)
  
-     ! multiplication factor according to slip system
-     if (lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID ) then
-       prm%omega    = prm%omega * 12.0_pReal      
-     else
-       prm%omega    = prm%omega * 8.0_pReal
-     endif
+     prm%omega                = config%getFloat('omega',  defaultVal = 1000.0_pReal) &
+                              * merge(12.0_pReal, &
+                                      8.0_pReal, &
+                                      lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID)
 
      ! expand: family => system
      prm%rho_mob_0    = math_expand(prm%rho_mob_0,   prm%N_sl)
@@ -808,8 +805,7 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
      if (dEq0(rho_dip_distance-rho_dip_distance_min(i))) then
        dot_rho_dip_climb(i) = 0.0_pReal
      else
-!       sigma_cl = norm2(matmul(Mp,prm%n0_sl(1:3,i))) ! ToDo: MD: correct?
-       sigma_cl = DOT_PRODUCT(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
+       sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
        if (prm%climbSwitch /= 0) then            
          b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i))
        else

From d2affaa414cab4169fc325ba2a42e18674902799 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 16 Aug 2019 14:04:12 -0700
Subject: [PATCH 2/2] clearer logic

---
 src/plastic_dislotwin.f90 | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90
index f1c42b678..9880ad868 100644
--- a/src/plastic_dislotwin.f90
+++ b/src/plastic_dislotwin.f90
@@ -118,7 +118,7 @@ module plastic_dislotwin
      sum_N_sl, &                                                                                    !< total number of active slip system
      sum_N_tw, &                                                                                    !< total number of active twin system
      sum_N_tr, &                                                                                    !< total number of active transformation system 
-     climbSwitch                                                                                    !< Switch to decide climb
+     climbSwitch                                                                                    !< Switch to decide climb !ToDo: Needs better name and should be a logical
    integer,                      dimension(:),     allocatable :: & 
      N_sl, &                                                                                        !< number of active slip systems for each family
      N_tw, &                                                                                        !< number of active twin systems for each family
@@ -280,7 +280,7 @@ subroutine plastic_dislotwin_init
      prm%SFE_0K               = config%getFloat('sfe_0k', defaultVal = 0.0_pReal)
      prm%dSFE_dT              = config%getFloat('dsfe_dt',defaultVal = 0.0_pReal)
      
-     prm%climbSwitch          = config%getInt('climbswitch',defaultVal = 1)
+     prm%climbSwitch          = config%getInt('climbswitch',defaultVal = 1) ! ToDo: turn on with config%key_exists('/name/')
  
      prm%omega                = config%getFloat('omega',  defaultVal = 1000.0_pReal) &
                               * merge(12.0_pReal, &
@@ -802,20 +802,22 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
        dot_rho_dip_formation(i) = 0.0_pReal
      endif
 
-     if (dEq0(rho_dip_distance-rho_dip_distance_min(i))) then
+     if (dEq0(rho_dip_distance-rho_dip_distance_min(i))) then !ToDo: use or here for climb switch (second if else not needed)
        dot_rho_dip_climb(i) = 0.0_pReal
      else
        sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
        if (prm%climbSwitch /= 0) then            
          b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i))
+         
+         v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) &
+              * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
+              
+         dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) &
+                              / (rho_dip_distance-rho_dip_distance_min(i))
        else
-         b_d = 0.0_pReal      
+         dot_rho_dip_climb(i) = 0.0_pReal
        endif
-
-       v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) &
-            * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
-       dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) &
-                            / (rho_dip_distance-rho_dip_distance_min(i))
+       
      endif
    endif significantSlipStress
  enddo slipState